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ABSTRACT 

We study the formation of H 2 in the ISM, using a modified version of 
the astrophysical magnetohydrodynamical code ZEUS-MP that includes a non- 
equilibrium treatment of the formation and destruction of H2. We examine two 
different approximations to treat the shielding of H 2 against photodissociation: a 
local approximation, which gives us a solid lower bound on the amount of shield- 
ing, and a method based on ray-tracing that is considerably more accurate in 
some circumstances but that produces results that are harder to clearly inter- 
pret. In both cases, the computational cost of determining H 2 photodissociation 
rates is reduced by enough to make three-dimensional high-resolution simulations 
of cloud formation feasible with modest computational resources. Our modifica- 
tion to ZEUS-MP also includes a detailed treatment of the thermal behaviour of 
the gas. 

In this paper, we focus on the problem of molecular cloud formation in grav- 
itationally unstable, initially static gas. (In a subsequent paper, we consider 
turbulent flow). We show that in these conditions, and for initial densities con- 
sistent with those observed in the cold, neutral atomic phase of the interstellar 
medium, H2 formation occurs on a timescale t > lOMyr, comparable to or longer 
than the gravitational free-fall timescale of the cloud. We also show that the 
collapsing gas very quickly reaches thermal equilibrium and that the equation of 
state of the thermal equilibrium gas is generally softer than isothermal. 

Finally, we demonstrate that although these results show little sensitivity to 
variations in most of our simulation parameters, they are highly sensitive to the 
assumed initial density n- x . Reducing ri\ significantly increases the cloud formation 
timescale and decreases the amount of hydrogen ultimately converted to H 2 . 
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Subject headings: astrochemistry — molecular processes — ISM: molecules - 
ISM: clouds 



1. Introduction 

Since essentially all observed Galactic star formation occurs within dense, self-gravitating 
molecular clouds, developing an understanding of the origin of these clouds is an important 
goal of research into star formation. Research in this area has typically focused on trying to 
answer a few basic questions: 

(i) How does the molecular gas form? In other words, what are the chemical processes 
involved? 

(ii) Where does the molecular gas form? Does it form before or after the assembly of the 
gas into dense clouds? 

(iii) How quickly does the molecular gas form? 



By far the largest constituent of the molecular gas is molecular hydrogen, H2, with other 
molecules such as CO being present only in small amounts, and so in practice the study of 
the formation of molecular gas is usually simply the study of the formation of H 2 . 

Although the chemistry of H2 formation in space remains an active field of study, the 
basic principles have been understood for some time. Gas-phase formation of H2 by direct 
radiative association is highly forbidden and proceeds at a negligible rate, while gas-phase 
formation via intermediat e molecular i ons such as H~ or H^ is strongly suppressed by the 
interstellar radiation field (lGloverll2003l ) and in any case cannot produce molecular fractions 
much higher than x-r 2 ~ 10~ 3 . Consequently, most Galactic H2 cannot have formed in 
the gas phase, leaving gra i n-sur face r eactions as the only v i able option. The pioneering 
work of iGould fc Salpeterl (119631 ) and iHollenbach fc Salpeterl (Il970l . Il97ll ) showed that H 2 
molecules could form on the surface of idealized dust grains with an effective rate coefficient 
compatible with that inferre d from UV observations of H 2 in the local interstellar medium, 



i?H 2 ~ 10 17 s 1 (|Juralll974l ). This remains widely accepted, at least for cold dus t, although 



there is ongoing debate about the efficacy of H 2 formation on warm dust (see e.g. iKatz et al. 
19991 : ICazaux & Tielensl[2004h . 



Answers to the other questions remain far more uncertain. One school of thought argues 
that H 2 forms in situ, in the locations presently occupied by the observed clouds. According 
to this picture, gas accumulates due to the action of large scale flows, which may be driven 



- 3- 



by large-scale grav itationa l instability (IKennicuttlll989l), magnetic instabi lities such as the 
Parker instability (jParkerl 1966, although see lKim. Ostriker fc Stond 120021 for a recent view 
of its importance), or may simply be part of the general turbulent velocity field, which 
itself is probably driven primarily by some combination of energy inpu t from supernovae 
and from the magnetorotational instability (IMac Low fc Klessenl 120041 ) . However, others 
have argued that the H 2 forms long before the molecular clouds themselves are assembled, 
residing in the interstellar medium in a diffuse st ate or in the for m of small cloudlets that 
eventually coalesce to form ob servable clouds (see lElmegreenlll990l and references therein, or 
Pringle. Allen fc Lubowll200ll for a more recent version of this model). Since coalescence will 
happen at a much faster rate in regions of converging flow, such as spiral arms, this model 
can be used to explain the enhanced star formation rates found within spiral arms. 

One way to discriminate between these explanations is by determining the ages of ob- 
served molecular clouds. If most molecular clouds are young, with ages comparable to their 
dynamical timescales, then this suggests that they are transient objects, and argues for a 
dynamical origin. On the other hand, if clouds are old, with lifetimes that are significantly 
greater than their dynamical timescale, then this is much easier to explain within a model 
in which clouds build up slowly and are dispersed slowly. 

In recent years, evidence that clouds are young has be en accumulating. For instance, 
Ballesteros-Paredes. Hartmann fc Vazquez- Semadenil (119991 ) argue that the absence of post- 
T Tauri stars with ages greater than 3 Myr in the Taurus-Auriga molecular cloud com- 
plex implies that the age of this cloud c o mplex can be no more than a few million years. 
Hartmann. Ballesteros-Paredes fc Berginl (120011 ) elaborate on this idea and show that simi- 
larly young ages are implied for most local star-forming regions. Additionally, the age dis- 
persion of stars in open clusters suggests that the molecular cloud complexes that give rise 
to them must have lifetimes of th e order of their d ynamical timescales, which are typically 
no more than a few million years (jElmegreenl 120001 ) . Short cloud lifetimes also make it eas- 
ier to understand the presence of supersonic magnetohydrodynamical turbulence within the 
molecu lar gas: in the absence of forcing, this will decay a way within a few turbulent crossing 
times ( jStone. Ostriker fc Gammid Il998l ; IMac Low! Il999l ). and so its presence in long-lived 
clouds requires there to be some form of external or internal driving, whereas its presence in 
short-lived clouds does not. 



Another key piec e of evidence for youthful clouds is discussed by iFukui et all (Il999l ) 
and lBlitz et all (120061 ). who show that giant molecular clouds in the Large Magellanic Cloud 
(LMC) are well-correlated with young stellar clusters, but do not correlate well with older 
clusters or with supernova remnants. They find that these very large clouds can only last 
~ 6 Myr before the onset of OB star formation, although they may last another 20 Myr, sup- 
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porte d by internal H II region expansion (e.g. lMatznerll2002l; iKrumholz. Matzner &: McKee 



Joung fc Mac Lowll2006l ). before 



20061 ) or external driving by background supernovae (e.g. 
their final dissolution by the same agents. 

However, for a model involving rapid cloud formation to be viable, it must be possible 
to produce the required quantity of H2 within a few million years. Given the relatively slow 
rate at which H 2 forms, it is natural to ask whether it is possible to satisfy this requirement. 
Simple back-of-the-envelope estimates made using the H 2 formation rate quoted above sug- 
gest that it is pos sible, provided that the mean density of the ma terial making up the cloud 
exceeds 10 3 cm -3 (IHartmann. Ballesteros-Paredes fc Bergmll200fl ). but we would ideally like 
to be able to confirm this result with more detailed num erical modeling. Various eff o rts in 



this direction have been made by a number of gr oups (jHennebelle &: Peraultl Il999l . 12000 



Koyama fc Inutsukall2000l . 120021 ; iBergin et a/.ll2004j ). but to date this modeling has generally 
been restricted to one or two dimensions, and has assumed an initially ordered, large-scale 
velocity field, such as a convergent flow, despite the observational evidence that the velocity 
field of t he neutral interstellar mediu m (ISM) is dominated by disordered, turbulent motions 



^see e.g. 



Lazarian Pogosyanl l2000l ) 



On the other hand, existing three-dimensional simulations of the neutral ISM, which do 



mal balance of the gas (e.g. 


KorDi et aiJl999: 


de Avillez 


2000; 


Wads 


2001 


Balsara et al. 


2004 




de Avillez & Breitschwerdt 


2004 




Slvz et al. 


2005: 



Kritsuk fc Normanl 2002a.b. 2004 



Joung &: Mac Low 



20061 ) have not previously included sufficient chemistry to follow the formation of H 2 and so 
have been unable to directly address the questions posed here. 

To bridge this gap, we have performed simulations of the neutral ISM using a hydro- 
dynamical code capable of following both the thermal balance of the gas and the formation 
and destruction of H 2 within it. Our goal is to determine how quickly significant quantities 
of H 2 can form in the dense neutral ISM, and to use this information to assess whether 
models in which cloud formation is assumed to be rapid are likely to work in practice. In 
this paper, we discuss in detail the numerical method used to follow the coupled chemical, 
thermal and dynamical evolution of the gas, and the tests used to verify the correctness 
of our implementation. We also present results from an application of our method to the 
problem of H 2 for mation in gas that is gra vitationally unstable, but not turbulent. In a 
companion paper (IGlover fc Mac Lowll2006l ; hereafter, paper II) we present results from a 
large suite of simulations that include the effects of supersonic turbulence. Although highly 
simplified, and probably not representative of real clouds, the non-turbulent models exam- 
ined in the second half of the present paper allow us to place an upper limit on the time 
required to form a molecular cloud, given neutral atomic gas with the density assumed in 
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our initial conditions. By comparing the results of these simulations with the results from 
the turbulent models examined in paper II, we can more clearly identify the effects of the 
turbulence, allowing us to demonstrate that supersonic turbulence significantly reduces the 
time required to form large quantities of H 2 . 

The structure of the current paper is as follows. In §[2] we describe the methods used to 
solve the equations of fluid flow, with a focus on our treatment of the thermal and chemical 
evolution of the gas, and in § [3] we present the results of some basic tests of our approach. 
In § m we describe and motivate the initial conditions used for our simulations. In § 
we present results from simulations of gravitationally unstable gas which is initially at rest, 
paying particular attention to the rate of H2 formation and the spatial distribution of the 
resulting molecular gas. Finally, in §E]we summarize our main results. 



2. Numerical method 
2.1. Magnetohydrodynamical equations 

The gover ning equations for the flo w of an inviscid, magnetized, self-gravitating gas can 



be written as (jStone fc Norman! Il992bl ): 



j£ = -,V.„, (!) 
= -Vj> - pV* + -^(V x B) x B, (2) 

an 

— = V x (v x B), (4) 

V 2 $ = AnGp, (5) 

where p e, p, v, B and $ are, respectively, the mass density, internal energy density, pres- 
sure, velocity, magnetic field and gravitational potential of the gas, where D/Di denotes the 
comoving derivative 

and where A denotes the net rate at which the gas gains or loses internal energy due to 
radiative and chemical heating and coolingQ Additionally, in a chemically reactive flow, 



Note that A > corresponds to a net loss of energy and A < to a net gain. 
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each chemical species satisfies an equation of the form 



Dp, 
Dt 



-PiV ■ v + Q - Dj 



(7) 



where pi is the mass density of species i, and where and D, represent its creation and 
destruction by chemical reactions. Finally, to close the set of equations, it is necessary to 
specify an equation of state relating the internal energy and the pressure. For an ideal gas 
we can write this as 

V = (7 - l)e, (8) 

where 7, the ratio of specific heats, depends upon the composition of the gas. For a gas with 
a number density n of hydrogen nuclei, a number density = O.ln of helium nuclei, and 
with a molecular hydrogen abundance xn 2 = In^jn and an electron abundance x e = n c /n, 
we can write 7 as 

3.3 + 3x c — 0.5xh 2 ' 

where we have assumed that the rotational degrees of freedom of H 2 are populated, but that 
the vibrational degrees of freedom are not. With the definition of xn 2 that we have adopted 
here, a value of xu 2 = 1.0 corresponds to gas in which all of the hydrogen is in molecular 
form, in which case 7 = 10/7. (Note that the more familiar 7 = 7/5 is for a gas which is pure 
H 2 ; the presence of the helium in our simulations causes a slight hardening of the equation 
of state). 

To solve this set of equations, we used a modified version of the publicly available ZEUS- 
MP hydrodynamical code. ZEUS-MP is a multi-phy sics, massively-parallel, message-passing 
code for astrophysical fluid dynamics (jNormanll2000l ). developed by the Laboratory for Com- 
putational Astrophysics at UC San Diego, which solves the equations of self-gravitating mag- 
netohydrodynamics (MHD) in three dim ensions. The algor ithms used in the ZEUS family 
of hyd ro codes are described in detail in IStone fc Norman! (19 92a,b) and lHawlev fe Stone 



Norman 



(1995). Their implementation within ZEUS-MP is discussed in iFiedlerl ( 119971 ) and 
(200 0J). The Poisson so lver used is a Fourier space solver that utilizes the FFTW library 
( iFrigo fc Johnson! Il998l ). Our modified version of ZEUS-MP is de rived from version 1.0b of 
the code. (For details of the more recently released version 2, see lHayes et a/.ll2006l ). 



We have modified ZEUS-MP in two main respects. Firstly, in order to follow non- 
equilibrium chemistry within the gas it is necessary to add an extra field variable for each 
chemical species that we wish to track. A natural choice of variable is the mass density 
of each species, as in that case in a medium with N non-equilibrium chemical species, we 
will have N equations of the form of equation [7] to solve. As discussed below in 
follow only two non-equilibrium species in the simulations presented here 



2~2l we 

H 2 and H+ - 
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and so have two such equations to solve. To solve these equations, we use operator splitting 
to separate the effects of advection (which is trea ted in the same fashion as advection of 



the total mass density; see lStone fc Normanlll992al ) from those of the chemical creation and 



destruction terms. In other words, during the reaction step we solve the equations 



dt 
dt 



Cn 2 (pn 2} pu+ } T) - D H2 (p H2 , Ph+, T), (10) 
C h +(Ph 2 ,Ph+,^) - Dh+(Ph 2 ,Ph+,^)- (11) 



The method of solution that we adopt is implicit finite differencing: we approximate equa- 
tions HOI & EH as 

Ph/ = Pk + ^ [Cg 1 - Dg 1 ] , (12) 

= Ph+ + At [C#+ — DjJ. 1 ] , (13) 

where the superscripts indicate values at the beginning and end of the timestep. The ad- 
vantage of using a first-order implicit method is that we can guarantee that the abundances 
will remain non-negative and that the solution will remain stable regardless of the size of 
the timestep chosen (although the requirement that we solve equations [T2] & [13] accurately 
still places some limits on the size of the timestep). The disadvantage of using an implicit 
method is that the resulting finite difference equations are coupled and must be solved it- 
eratively. Moreover, the fact that the creation and destruction terms also depend strongly 
on the internal energy of the gas (through the temperature T) suggests that we should solve 
these equations simultaneously with the energy equation. We therefore defer discussion of 
the solution of equations [T2] & [T3] until after we have discussed the modifications that we 
have made to the treatment of the internal energy equation. 

This has been modified to include a term representing the combined effects of radiative 
and chemical heating or cooling, i.e. the A term of equation [31 Details of the processes 
included are given in § 12.31 below and are summarized in Table [2J Solution of the resulting 
equation proceeds much as it does in the unmodified version of ZEUS-MP: the equation is 
operator split, with the effects of artificial viscosity, compressional heating and advection 
treated separately. Our treatment of the artificial viscosity an d advection steps mirrors that 



in the unmodified version of the code, as discussed in detail in lStone fc Normanl (Il992al ); we 
will not discuss it here. However, during the compressional heating step, instead of solving 
the equation 

de 

— = -pV ■ v, (14) 

we solve 

3c 

— = - p v • v - A. (15) 
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To solve equation [T5], we use an algorithm based on a combination of the ZEUS-3D pdv 
cool algorithm (originally implemented by M. Norman and subsequently m odified by M.- 
M. M ac Low, J. Stone and D. C larke) with the i mplici t algo rithm used by ISuttner et al. 
r ll997h and further developed by IPavlovski etall J2002h and ISmith fc Rosenl J2003h . We 



construct the following implicit approximation to equation 



e i+i _ e i 
At 



-pV ■ v l - A m , 



(16) 



where p represents the time-centered pressure, which we approximate as p ~ 0.5[p(e 4+1 ) + 
p(e 1 )], and where A l+1 is the cooling rate at the end of the timestep. We rearrange this 
equation to give 

;i-g)e'-A i+1 At 



1 + q 



(17) 



where 



At 



? = T (7"1)(V-^), 



and then solve equation [T7] together with equations [12] & [13] using a form of Gauss-Seidel 
iteration. This works as follows: 



1. Update pn 2 ) using old values for p H + and e. 

2. Update pn+, using new value for pn 2 but old value for e. 

3. Update e, using new values for p H2 an d Ph+- 

4. Test for convergence. If not converged, return to step 1. 



On our first pass through the loop, we take the old values for p H + and e to be those at the 
start of the current timestep. On subsequent passes through the loop, we instead use the 
values from the previous iteration. 

Since the internal energy converges more slowly than either of the chemical abundances, 
we monitor its value to determine when to stop the iteration, halting once the relative 
difference between updated and old values becomes less than 1CU 7 . On rare occasions, the 
iteration may fail to converge. In this case, we switch to solving for e using a more expensive 
bisection algorithm that is sure to find a solution. In this case, we lose the benefits of the 
iteration for refining our initial guesses for pn 2 and pn+- Fortunately, the abundances of 
H2 and H + generally do not vary much during a single timestep and so the error that this 
introduces is probably not significant, particularly since the iteration almost always succeeds. 
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To co mpute the cooling timestep, At coo i, we use the same procedure as in lSuttner et al. 



( 119971 ) and iPavlovski et all (120021 ): we compute A and then require that the timestep be 



such that the internal energy will not vary by more than 30% given this value of A, i.e. 

Atcooi = 0.3^. (19) 

To improve the efficiency of the code, we also make use of subcycling. Rather than 
constraining the global timestep of the simulation At to be less than At coo i, we instead use 
the same global timestep as we would within the unmodified version of the code, but solve 
the chemistry and cooling substeps alone over the shorter cooling timescale, repeating the 
procedure for as many times as is necessary to make the total elapsed time equal to At. 
Thus, if Atcooi > At, we proceed through the chemistry and cooling substep only once per 
simulation timestep, while if At coo i <C At, we do so many times, using the updated values 
of the chemical abundances and internal energy from the end of one substep as the input 
to the next. We subcycle at the level of the individual grid zones, so only zones for which 
Atcooi < At take multiple chemistry and cooling substeps per hydrodynamical timestep. The 
gain in computational efficiency from this subcycling procedure is difficult to quantify, as it 
depends on the physical state of the gas, but tests indicate that it can easily be as much as 
an order of magnitude. 



2.2. Chemistry 



The chemical composition of the IS M is complex. Over 120 different molecular species 
have been detected in interstellar space (jWiklindll2004l ) and while many of these are found 
in detectable amounts only in dense, well-shielded gas, there rem ain a significant number 
that have been detected in diffuse, unshielded gas (see, for instance lO'Neill. Viti fc Williams 
20021 . and references therein). A full chemical model of the ISM can easily involve several 
hundred different atomic and molecular species a nd several thousand different reactio ns, even 
if reactions on grain surfaces are neglected (e.g. 



Le Teuff. Millar fc Markwickll2000h 



It is currently impractical to incorporate this amount of chemistry into a 3D hydro- 
dynamical code such as ZEUS-MP, due to the extreme impact it would have on the code's 
performance. Fortunately, much of this chemistry is not relevant to our current study, and 
can be neglected without significantly compromising our results. Recall that the main goals 
of this paper are to determine the timescale on which molecular clouds form in gravitation- 
ally unstable atomic gas, and the fraction of this gas that is converted to molecular form. 
To achieve these goals, we clearly need to be able to follow the formation and destruction 
of H 2 with a reasonable degree of accuracy. Beyond this, however, the only chemistry that 
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we really need to be concerned with is that which plays a role in determining the thermal 
balance of the gas. In other words, we need only follow the chemistry of H 2 , and of a few 
other major coolants such as C + or O; the chemistry of other molecules such as CH, while 
undeniably of interest, is not central to the task at hand and so can be sacrificed in pursuit 
of efficiency. 



T o iden tify the major coolants in the diffuse ISM, we can use the results of lWolfire et al. 



(119951 . 120031 ). who studied its thermal evolution in some detail. They show that in the warm 
neutral medium (WNM), which has a characteristic temperature T ~ 8000 K, most of the 
cooling comes from Lyman-a emission from atomic hydrogen, electron recombination with 
small grains and polycyclic aromatic hydrocarbons (PAHs), and fine structure emission from 
atomic oxygen. Cooling in the cold neutral medium (CNM), which has a characteristic 
temperature T < 300 K, is dominated by fine structure emission from ionized carbon, C II , 
with fine structure emission from oxygen also contributing significantly in the warmer parts 
of the CNM. These two coolants also dominate in the thermally unstable temperature regime 
300 < T < 8000 K. Heating in all three regimes is dominated by photoelectric emission from 
dust grains and PAHs; the heating rate is therefore primarily determined by the strength 
of the ultraviolet background and by the dust to gas ratio, but it also has a dependence on 
the electron abundance (see equation H3] in § 12.31) . We also include fine structure cooling 
from Si + . Although never the dominant coolant, Si + does produce ~ 10% of the total fine 
structure emission at temperatures T > 200 K at all densities, and can produce as much as 
30% of the emission at high densities, owing to the relative ly large size of Einstein coefficient 



for the Si+ 2 P 3/2 -> 2 P 1/2 transition flSilva fc Viegasll2002r ) 



In denser gas that is shielded from the ultra violet background, other coo l ants such as CO 



and H ?0 be come important, as demon s trated by lGoldsmith fc Langerl (119781 ) . iNeufeld fc Kaufman 



(119931 ) and INeufeld. Lepp fc Melnicki (119951 ). However, since we are primarily interested in 
the transition from atomic to molecular gas, and since these coolants will not become im- 
portant until after the gas is already fully molecular, we have chosen to make another major 
simplification in the chemistry: we assume that carbon, oxygen and silicon remain primarily 
in the form of C II , O I and Si II respectively, throughout the simulation. By making this 
assumption, we essentially reduce the chemistry that must be followed to that of only three 
species - electrons, neutral hydrogen and H 2 - at the cost of computing an incorrect cooling 
rate in dense molecular gas. The dense gas in our simulations is also typically rather cold, 
with a temperature of no more than a few tens of Kelvin, and so we will probably under- 
estimate the cooling rate within it, since in these conditions molecular coolants such as CO 
are more effective than the C II , O I and Si II included in our simulations. This means that 
we probably overestimate the temperature of this gas, but the fact that the dense gas is so 
cold shows that we do not overestimate its temperature by a large amount, and we therefore 
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do not expect this approximation to have a significant impact on the H2 formation rate. 



Having discussed the underlying assumptions, we now present our chemical model. We 
adopt standard solar abundanc es of hydrogen and hel ium, along with abundances of carbon, 



oxygen and silicon taken from ISembach et al\ ( 120001 ) : xq = 1.4 x 10 



Xq 



3.2 x 10' 



xsi = 1.5 x 10~ 5 . We assume that helium remains neutral and plays no direct role in the 
chemistry, and that carbon and silicon remain singly ionized throughout the simulation. The 
ionization state of oxygen is assumed to track that of hydrogen, due to the influence of the 
charge transfer reactions 



0+ + H -> + H+, 
H+ + -> H + 0+ 



(20) 
(21) 



As previously noted in § 12.11 we follow directly the abundances of two chemical species, 
namely H + and H 2 . The abundances of the other major species - atomic hydrogen and 
electrons - are computed from the conservation laws: 

xr = x H ,tot - %h+ ~ x H2 , (22) 

and 

x e = x H + + x c + + x S i+, (23) 

where x# itot is the total abundance of hydrogen nuclei, in all forms, and where x c + and x Si + 
represent the abundances of ionized carbon and silicon respectively. We neglect the small 
contributions to x e made by electrons from other ionized metals, such as Mg or S, since these 
are small compared to the contribution from carbon. We also neglect any contribution from 
+ , as this will be negligible compared to the contribution from hydrogen. 

Our assumption that carbon and silicon remain ionized throughout the simulation be- 
comes inaccurate in dense, well-shielded gas, where UV photoionization of neutral carbon 
and silicon becomes ineffective. We therefore overestimate the fractional ionization of gas 
in these regions. However, our main motivation for tracking the fractional ionization is to 
compute the photoelectric heating rate accurately, and since this is unimportant in dense, 
shielded gas, the inaccuracy in the fractional ionization is likely unimportant there. 

The reactions included in our chemical model are summarized in Table [TJ In most cases, 
we also list the source or sources from which we took the rate coefficient used in our model. 
The exception is H 2 photodissociation, the treatment of which is discussed in detail in § 12.2.11 
below. Two other reactions deserve further comment: the collisional dissociation of H 2 by 
atomic hydrogen 

H 2 + H -> 3H, (24) 
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and molecular hydrogen 



H 2 + H 2 — > H 2 + 2H. 



(25) 



The reaction coefficients for both of these reactions are density-dependent, since they are 
sensitive to the population of the vibrational and rotational levels of H 2 . To treat the former, 
we use a rate coefficient 



log k H 



n/n c 



1 + n/n c 



log few] 



1 + n/n c 



log few; 



(26) 



where fc^i is the l ow de nsity limit of the collisional dissociation rate and is taken from 
Mac Low &: Shulll (119861 ). while fcn,h is the high density limit, taken from iLepp fc Shull 
( 19831 ) . The critical density, n cr , is given by 



n r 



+ 



XR 2 



n. 



cr,H 2 



(27) 



where n CTj H and n CTj H 2 are the critical densities in pure atomic gas with an infinitesimally 
dilut e quantity of H 2 and in pure molecular gas respectively. The first of these values is taken 
fro mlLepp fc Shulll (I1983h. but has be en decreased by an order of ma gnitude, as recommende d 
by iMartin. Schwarz fc Mandyl (119961 ) ; the other value comes from IShapiro fc Kangi (119871 ) . 
To treat the collisional dissociation of H 2 by H 2 we use a similar expression 



log k 



H 2 



n/n c 



1 + n/n c 



log A; H2 ,h + 



1 + n/n c 



log k 



H 2 ,h 



(28) 



where the low density limit, A;h 2 ,i, is t aken from IMartin. Keogh fc Mandyl (119981 ) and the 
high density limit, ku 2 ,h, is taken from lShapiro fc Kangi (119871 ). The collisional dissociation 
rates computed in this way are acceptably accurate when wh nn 2 or <C nu 2 , but may 
be less accurate in gas with n H ~ nn 2 ; further study of the collisional dissociation of H 2 in 
gas which is a mixture of H and H 2 would be desirable to help remedy this. 

We do not include the gas phase formation of H 2 via the H~ or ions, as in the typical 
conditi ons of the Ga lactic ISM, this is unimportant compared to the formation of H 2 on dust 
grains (jGloverl 120031 ). 



As can be seen from Table (U our H + chemistry is straightforward. However, a couple 
of reactions are deserving of comment. One is gas-phase recombination, where we note that 
we use the case B value for the recombination coefficient, as this is the most appropriate 
choice in all but the most highly ionized gas. The other is the cosmic ray ionization of 
H. In the majority of our simulations, we use a value of the cosmic ray ionization rate 
( = 10~ 17 s" 1 . This is consistent (within the error b ars) with recent de t erminations of the 



ionization rate in dense cores within molecular clouds (jCaselli et a/.lll998l ; iBergin et a/.lll999 
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van der Tak van Dishoeckll2000l ). However, measurements of the Hjj" column density along 



sight-lines passing through diffuse gas are only explained by a much larger value of C, within 



the r ange 1(T 16 < C < 10~ 15 s" 1 flMcCall et al.ll2003l : iLisztl 120031 : iLe Petit. Roueff fc Herbst 



2004 ). where the uncertainty is due primarily to the uncertain distribution of gas along 



the observed lines of sight. The cause of this discrepancy is currently unknown (although 



Padoan fc Scald 120051 suggest one possible mechanism), and so for simplicity we take ( to 
be constant, independent of the gas density. We briefly examine in § 15.41 the consequences 
of adopting a larger value for £. 



2.2.1. H2 photodissociation 



Following iDraine &: Bertoldil (119961 ). we can write the photodissociation rate of H2 in 
optically thin gas as 

k phfi = 3.3 x 10- n x s~\ (29) 



where we have assumed that the ultraviolet field has the same spectral shape as the IDraine 
( 119781 ) field, and where x is a dimensionless factor which characterizes the i ntensity of the 
field at 1000 A relative to the iHabingl ( 19681 ) field; note that for the original [DraineJ ( 19781 ) 
field, X = 1-7. 

By balancing this dissociation rate against the rate at which H 2 forms on dust grains, 
we can easily show that the equilibrium H 2 fraction in cold gas is given by 

T \ 1/2 



Xh 2 



10 



"V 1 



100 



(30) 



which is clearly far less than one unless % is very large. Since observatio ns indicate that 
signifi cant molecular gas is present at densities below % 2 < 10 4 cm~ 3 (see e.g. lFalgarone et al. 



19981 ). it is clear that some shielding of interstellar H2 molecules from the effects of UV 
photodissociation must occur, and that simulations of molecular cloud formation which do 
not take this shielding into account are going to be tremendously inaccurate. 

H2 can be shielded against photodissociation in two main ways: by line absorption due 
to other H 2 molecules (i.e. self-shielding), and by continuo us absorption due to dus t. The 
effects of both processes have been treated in some detail by IDraine &: Bertoldil (119961 ). They 
show that in shielded gas the photodissociation rate can be written approximately as 

k ph = /shieid(^ H2 )e- Td - 1000 A; ph ,o, (31) 

where /shield is a numerical factor accounting for the effects of self-shielding, r^iooo is the 
optical depth due to dust at IOOOA, and fc p h,o is the unshielded photodissociation rate, given 
by equation [29] above. 
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Draine Bertoldil show that for a static, plane parallel slab of gas, /shield is well approx- 



imated by 

0.965 0.035 r ,,,,, , , 

/shield = ITT1 ^ + exp [-8.5 x 10-(1 + X )V] , (32) 

where x = N-g_ 3 /5 x 10 14 cm~ 2 with Nn 2 being the H 2 column density, and b 5 = fe/10 5 cms -1 , 
where b is the Doppler broadening parameter. 

If the gas is not static or in uniform motion, but instead has a spatially varying velocity, 
then equation [32] will overestimate the amount of self-shielding that occurs. The reason 
for this is that if the velocity field is non-uniform, then the relative velocity between any 
two fluid elements will, in general, be non-zero. Consequently, the contribution to the total 
absorption coming from the first fluid element will be Doppler shifted when viewed from 
the rest frame of the second fluid element. If this Doppler shift is large compared to the 
line widths of the Lyman- Werner band transitions, then the effect is to significantly reduce 
the extent to which the absorption coming from the first fluid element contributes to the 
self-shielding seen by the second fluid element. For H2 column densities Nr 2 < 10 17 cm" 2 , 
the intrinsic widths of even the strongest Lyman- Werner band transitions are unimportant 
and the line profiles are dominated by Doppler broadening. In this regime, the neglect of 
changes in the velocity along the line of sight is justified if the differences in velocity are 
much smaller than the sound speed of the gas, and equation [321 remains a good estimator for 
/shield- On the other hand, if the differences in the velocity are much greater than the sound 
speed, as will generally be the case in a supersonically turbulent flow, then equation [32] will 
significantly overestimate the degree of self-shielding that will actually occur. 

For A^h 2 > 10 17 cm -2 , the line widths of the strongest transitions are dominated by 
Lorentz broadening, and the degree of self-shielding in these lines becomes insensitive to the 
velocity distribution of the gas, unless the range of velocities is extremely large. However, the 
total H2 photodissociation rate remains sensitive to the velocity dispersion, as a large fraction 
of the total rate is due to absorption in weaker lines, whose widths are still dominated by 
Doppler broadening. Only for Nr 2 > 10 18 cm -2 does the total photodissociation rate become 
relatively insensitive to the velocity distribution of the gas. 

An additional problem with using equation[32]to compute the self-shielding factor is that 
it requires us to compute the H2 column densities along a large number of lines of sight for 
every zone in our simulation. Unfortunately, the computational cost of doing so is very large: 
for a simulation with N fluid elements, the cost of solving for the column densities scales as 
0(iV 5 / 3 ) (assuming that we require an angular resolution that is well matched with the spatial 
resolution of the code), compared to scalings of O(N) for the hydrodynamical algorithms 
and 0(N log N) for the self-gravity. Consequently, solving for the column densities would 
quickly come to dominate the cost of the simulation, and would prevent us from performing 
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any high resolution simulations. Moreover, while parallelization of the code would help to 
some extent, efficient parallelization of the radiative transfer is difficult, owing to the way in 
which it couples together widely separated regions in the computational volume. 

Computation of the optical depth due to dust, r^iooo, is in principle considerably sim- 
pler than computation of /shield, since it involves continuous absorption rather than line 
absorption, an d so Doppler effects can generally be neglected (at least for the velocities of 



interest here). ICardelli. Clayton fc Mathisl (119891 ) showed that it is possible to approximate 



the extinction observed along many lines of sight in the Galaxy by a one-parameter family 
of curves, where the controlling parameter is Ry, defined as 

where Ab and Ay are the extinctions in the B and V bands respectively. For R y = 3.1, 



which is typical for many lines of sight in the diffuse ISM. lDraine fc Bertoldil (119961 ) quote an 
effective attenuation cross-section for the dust of o" d lo0 o = 2 x 1CT 21 cm 2 . The corresponding 
dust opacity is given by 

r dilo0 o = 2 x l(T 21 iVH, tot , (34) 

where iVfi^ot = + + 2A^h 2 is the total column density of hydrogen nuclei between 
the point of interest and the source of the radiation. Therefore, to compute the effects of 
the dust shielding, we again need to compute appropriate column densities, although in this 
case the effect of velocity differences along the line of sight is unimportant. As in the case of 
self-shielding, the cost of solving for the column densities is 0(N 5 ^ 3 ) for a simulation with 
N fluid elements, and so a high resolution treatment of both the fluid flow and the radiation 
field is not computationally feasible. 

To overcome these difficulties, we have performed simulations using two different approx- 
imations. Our first approximation is very simple. We assume that the dominant contribution 
to the shielding of a given fluid element comes from gas in the immediate vicinity of that 
element. In the case of H2 self-shielding, this assumption can be justified to some extent by 
the fact that the relative velocity between this gas and the point of interest will typically 
be small, whereas gas at larger distances will typically possess a much larger relative veloc- 
ity, particularly in supersonically turbulent gas. Additionally, in gravitationally collapsing 
regions, the local gas density will be substantially larger than the density in most regions of 
the flow, further increasing the importance of the spatially local contribution. In the case 
of dust shielding, the former justification no longer holds true, although the latter remains 
valid. To give a quantitative form to this approximation, we note that we continue to use 
equation [32] to compute the self-shielding factor, but use an H 2 column density given by 

Ax 

A% 2 = — n Ha , (35) 
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where n^ 2 is the H2 number density in the zone of interest and Ax is the width of the zone, 
measured parallel to one of the coordinate axes. Similarly, we continue to use equation [3U 
to compute the dust opacity, but replace N^ )tot with the local approximation: 

Ax 

Nn t tot = — {nn + n H + + 2n U2 ) , (36) 

where %, n H + and nn 2 are the local number densities of atomic hydrogen, H + and H 2 
respectively. In the remainder of this paper, we refer to this approximation as the local 
shielding approximation. 

The main advantage of the local shielding approximation is that we can be certain 
that we are underestimating the true amount of shielding and hence overestimating the 
H 2 photodissociation rate. Consequently, H 2 fractions computed in simulations which use 
the local shielding approximation will be strong lower limits on the true values, and the 
corresponding H 2 formation timescales will be solid upper limits. Moreover, in paper II we 
show that this very simple approximation proves to be surprisingly accurate in supersonically 
turbulent gas, although it is significantly less effective in the simulations of smooth collapse 
presented in this paper. 

A significant disadvantage of the local shielding approximation is that it makes the 
photodissociation rate dependent on Ax and hence on the resolution of the simulation. 
Consequently, the abundance of H 2 that we expect to find in gas in our simulations that is 
in photodissociation equilibrium also becomes resolution-dependent. To see this, consider 
the equation for the equilibrium H 2 abundance: 

^H 2 ,eq 2i?f orm 

-n, (37) 



1 3'H2,eq Rph 



where -Rf or m is the formation rate of H 2 on dust grain surfaces, taken from lHollenbach fc McKee 



(119791 ). and -Rph = /c p h^H 2 - Since _R pri depends on Ax but Rf OIia does not, the equilibrium 
H 2 abundance xu 2 ,eq must inevitably depend on Ax. How significant this is for any given 
grid zone depends in large part on the total hydrogen column density and the H 2 column 
density for that zone. If both are sufficiently small that the zone is optically thin to UV 
radiation, then changing the size of the zone will have little effect on R p ^, provided that 
the zone remains optically thin. Conversely, if either A^H,tot or Ah 2 f° r that zone are large 
enough to make the zone highly optically thick to H 2 photodissociation, then changing the 
zone size will have a large effect on Rph, but only a very small effect on XH 2 , e q, which will 
remain close to one. However, for zones that lie in the intermediate regime between these 
two limiting cases, the effect on XH 2 , eq can be significant and it is here that our results will 
be least accurate. 
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The second approximate technique that we have used to compute the photodissociation 
rate again makes use of equations [32] and [3H However, in this approach we compute the H 2 
and dust column densities seen by each grid zone by averaging over a small number of lines of 
sight. Specifically, we compute column densities in both the positive and negative directions 
along each of the three coordinate axes of the simulation, compute the associated values of 
/shield and Td,iooo f° r eac h of these lines of sight, and then compute the total photodissociation 
rate by means of a suitably weighted average of these values. The main advantage of this 
approach is that its computational cost should not be very much greater than the cost of 
solving the hydrodynamic equations, as in both cases the cost of the algorithm is O(N). 
Moreover, by restricting the lines of sight that are treated to be those parallel to the coor- 
dinate axes, we also limit the amount of communication between different processors that is 
required, and so limit the adverse impact on the scalability of the code. An approach of this 
type has previously been used to study the stab i lity a nd dynamics of low mass molecular 



clouds and Bok globules (INelson fc Langerl 119971 . Il999l ) , while a very similar approach has 



been used to stud y the formation of H 2 in early protogalaxies illuminated by an intergalactic 



UV background (jYoshida et al.1 120031 ). We refer to this approximation in the remainder of 



this paper as the six-ray shielding approximation. 

An obvious disadvantage of this approach is the extremely coarse angular resolution of 
the radiation field that it provides. This poor angular resolution will cause us to overesti- 
mate the amount of shielding in some regions, and underestimate it in others: the precise 
details will depend on the particular form of the density field, but in general we will tend 
to underestimate the amount of shielding whenever the volume filling factor of dense gas is 
small. Another significant problem is that this approach does not take account of velocity 
structure along any of the lines of sight. It therefore may significantly overestimate /shield 
in a supersonic flow, particularly if the gas is turbulent. For the main problem that we 
are interested in investigating - the determination of the H 2 formation rate in dynamically 
evolving, cold atomic gas - this is problematic, as it may lead us to derive an artificially short 
timescale for H2 formation. Nevertheless, in many scenarios (including the simulations of 
smooth collapse presented later in this paper), the six-ray shielding approximation produces 
significantly more accurate results than the local shielding approximation. 

Another possible approach to treating the effects of shielding which we considered but 
discarded makes use of the fact that the contribution towards the total H 2 photodissociation 
rate made by any individual Lyman- Werner line is directly proportional to the penetration 
probability for that line. For a line of sight n, this is given by 



Pp(n) = / 4>(u)e~ T ^du, (38) 
Jo 

where <j>{y) is the line profile function for the line in question and r(u, n) is the optical depth 
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at frequency v in the direction n. This penetration probability is analogous to the more 
commonly encountered escape probability, and this similarity can be fruitfully exploited. In 
particular, if the gas h as a large mon otonic velocity gradient in the direction n, then the 



Sobolev approximation (jSobolevlll957l ) can be used to compute -P p (n) given only the local H 2 
number density and the magnitude of the velocity gradient. If the velocity gradient is large 
in all directions, then we can repeat this procedure for many lines of sight, and with suitable 
averaging can derive a mean penetration probability for the line. Finally, by repeating this 
for each of the lines which contribute to the total photodissociation rate, we can compute 
the rate itself, using only local quantities. 

There are, however, two major drawbacks to this approach. First, it can only be relied 
upon to give accurate answers in conditions where the Sobolev approximation applies. This 
formally limits the applicability of this method to regions where the velocity gradient is 
monotonic and where any variations in density, temperature or H2 abundance occur far more 
slowly than variations in velocity. These conditions are not satisfied in turbulent molecular 
clouds, and so although the S obolev approximation can sometime be fruitfully applied (see, 



for instance lOssenkopf 119971 ). in general it will not give accurate results. Moreover, in 
contrast to our local shielding approximation, we cannot be confident that we know the 
sense of the inaccuracy, i.e. whether we produce too much or too little H 2 , as this will 
depend on the physical conditions within a given simulation, and may also change over time 
within that simulation. Therefore, although this approach avoids the resolution dependence 
of our local shielding approximation, it does so at the cost of producing results that have 
no clear interpretation, being neither lower nor upper limits. The other major drawback to 
this approach is that it cannot be used to model dust absorption, as in this case the Sobolev 
approximation simply does not apply. 

To sum up, our local shielding approximation has the advantages of simplicity and a 
straightforward interpretation, while our six ray shielding approximation will often be more 
accurate but has a less clear interpretation. Although neither approximation is ideal, by 
comparing and contrasting the results from both, we can draw important conclusions about 
the formation timescale of H 2 in static and in turbulent gas. 



2.3. Heating and cooling 

In order to solve equation [151 we need to calculate A, the net rate at which the gas 
gains or loses energy due to radiative and/or chemical heating and cooling. We can write A 
as the sum of a heating and a cooling term: 

A = A cool — r hcat . (39) 
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As discussed in § I2.2[ a number of different processes contribute to A coo i- At high temper- 
atures (T > 8000 K), much of the cooling is provided by excitation of the resonance lines 
of hydrogen, helium or heavier elements. To treat cooling in thi s temp erature range, we 
adopted a tabulated cooling function from ISutherland k, Dopital (119931 ): specifically, the 
cooling function listed in Table 10 of their paper, which assumes [Fe/H] = -0.5. Significant 
cooling also comes from the recombination of elec trons with small dust grains and PAHs. We 
incorporate thi s using a cooling rate t aken from lWolfire et all ( 120031 ). based on an original 
formulation by lBakes fc Tielend (119941 ): 



A r 



4.65 x 10" 3 V pah T a9 Wn 



ergs cm 3 s , 



where (3 = 0.74/T - 068 , ip is given by 



XesT 



1/2 



(40) 



(41) 



X represents the strength of the UV background in the gas after dust 



where Xefr = e" 

shielding is taken into a ccount (see iBergin et all 120041 ). and where pa h is an adjustable 
parameter introduced by I Wolfire et al\ (120031 ) to make the heating rate consistent with the 
values of the electron attachmen t and electron rec ombination rates that are inferred obser- 
vationally for PAHs; the original iBakes fc Tielend treatment corresponds to pa h = 1. 



At lower temperatures, the contribution of these coolants becomes negligible, and cool- 
ing by C II and O I fine structure lines domin ates. To compute the cooling from C II , we 
used atomic data from ISilva &: Viegasl (120021 ) , together wi th collisional de-excitation r ates 
from iFlower fc Launayl (119771 ) for collisions with H9 . from iHollenbach &: McKed (119891 ) for 
collisions with atomic hydrogen at T < 20 00 K, fromlKeenan et all Jl986h for collisions with 



atomic hydrogen at T > 2000 K, and from I Wilson fc Belli (120021 ) for collisions with electrons. 



For O I fine structure cooling, we used atomic data from lSilva &: Viegad (120021 ). together 
with collisional rates provided by D. Flower (priy . comm .) for collisions with H and H2, 
as well as rates from iBell. Berrington &: Thomad (119981 ) for collisions with electrons and 
Pequignotl (1990, 1996) for collisions with protons. 



In addition to cooling from C II and O I , we also included contributions from Si II fine 
structure line emission - which can be more effective than C II cooling if the temperature 
and density are both large - and from the rotational and vibrational lines of H2. 

To compute the Si 1 1 cooling rate, we again used atomic data from Silva fc Viegasl (2002h . 



and co llisional rates from lRouefil (119901 ) for collisions with atomic hydrogen and from lDufton fc Kingston 
(Il99ll ) for collisions with electrons. De-excitation rate coefficients for collisions between Si II 
and H 2 were unavailable, and so were arbitrarily set to zero; however, this is unlikely to 
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introduce a significant error into the computed cooling rate in our simulations as gas with a 
significant molecular fraction is typically far too cold for Si II cooling to be effective (see, for 
instance, § 15.41) . 



For cooling due to H 2 , we use the cooling function from lLe Bourlot. Pineau des Forets fc Flower 



( 119991 ) . which we have extended to temperatures below 100 K by assuming that only the 
J = 2 — > and J = 3 — > 1 transitions contribute significantly to the cooling rate. For sim- 
plicity, we also fix the ortho:para ratio at 3:1. However, variations in this ratio are unlikely to 
significantly affect the H 2 cooling rate at temperatu res at which it contributes significantly to 
the t otal cooling rate (see, for instance, Figure 5 in lLe Bourlot. Pineau des Forets & Flower 
1999h . 



In each case, we assumed that cooling occurs in the optically thin limit. This is a 
reasonable assumption in diffuse gas, but breaks down in dense, high column density cores. 
However, in these conditions of high density and high optical depth, we would expect all of 
the hydrogen to already be in molecular form (and indeed much of the carbon and oxygen 
to be in the form of CO, rather than Clland Olas assumed here) and so errors in the gas 
temperature in these dense cores are unlikely to have any significant effect on our results. 

In addition to these processes, we also include t he effects of collisional transf er of energy 
between gas and dust, following the prescription of iHollenbach fc McKed (119891 ) . This acts 
to cool the gas whenever T gas > T dust , and to heat it when T dust > T gas . However, this is not 
an important process at the gas densities studied here, although it does become increasingly 
important in higher density gas. 

Finally, we also allow for the effects of cooling due to the collisional dissociation of H 2 
and collisional ionization of H, although in practice neither process is of much importance 
in our simulations. 

As far as heating is concerned, the most important contribution to Theat comes from the 
heatin g produced by ph o toelec tric emission from UV-ir radiated dust grains and PAHs. Fol- 
lowing |Bakes_^IIlelenj (119941 ) and lWolfire et al\ (120031 ) , we write the photoelectric heating 
rate as 

r pe = 1.3 x 10" 2 Vxcfr ergs s~ x cm" 3 , (42) 
where e is the heating efficiency, given by 



4.9 x 10~ 2 



+ 



3.7 x 10~ 2 (T/10 4 ) ' 7 



i + 4.o x io- 3 ^ - 73 1 + 2.0 xio-y 

with ip as given by Equation [4X1 above. 

Additional contributions to I^cat come from several other processes: 



(43) 
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Photodissociation of H2 by far-ultraviolet (FUV) radiation 



Following iBlack fc Dalgarnol (Il977l ). we assume that each photodissociation deposits 
0.4 eV of heat into the gas, giving us a heating rate 



r ph = 6.4 x 10 13 k p hnn 2 ergs s 1 cm 



(44) 



where fc p h is given by Equation I3T1 above. 



Excitation of H 2 by FUV radiation 

As well as dissociating some of the H 2 , the FUV radiation also produces vibrationally 
excited H 2 via radiative pumping. At high densities, this leads to heating of the 
gas, as most of the excited molecules undergo collisional de-excitation. We adopt 
a radiative pumping rat e that is 8.5 times larger than the photodissociation rate 
(jDraine fc Bertoldil 1 19961 ). and assume that each excitation trans fers an average of 
2 (1 + ricr/n)" 1 eV to the gas (IBurton. Hollenbach fc Tielensl Il990l ). where n cr is the 
critical density at which collisional de-excitation of vibrationally excited H 2 occurs at 
the same rate as radiative de-excitation. As previously noted, our value for Ti cr is a 
weighted harmonic mean of the value for H 2 -H collisions given bvlLepp fc Shulll (119831 ) 



and the value for H 2 -H 2 collisions given by [Shapiro fc Kangl (119871 ) . 



H 2 formation on dust grains 

The formation of an H 2 molecule from two hydrogen atoms releases approximately 
4.5 eV of energy. In all likelihood, some fraction of this energy will go into the 
rotational and vibrational excitation of the newly-formed H 2 molecule, with the re- 
mainder going into the t ranslational energy of the molecule or into heating the grain 
(jDuley fc Williams! Il993l ). The fraction of the total energy going into rotational an d 
vibrational excitation remains a subject of investigation (see e.g. iRoser et q/.ll2003l ). 
but in our simulations we assume that this fraction is close to one. As with FUV 
pumping, most of this energy will be lost through radiative de-excitation of the H 2 
molecule if n <C n cr , but will be converted to heat if n 3> n CT . We therefore adopt an 
H 2 formation heating rate of 



r H2 = 7.2 x 10 



-12 



R 



H 2 



[1 +n cv /n) 



ergs s 1 cm 3 , 



(45) 



where Rn 2 is the rate of H 2 formation on dust grains and n cr is the critical density, 
calculated as described above. 



• Cosmic ray ionization 
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Following iGoldsmith fc Langerl (119781 ). we assume that each ionization deposits 20 eV 
of energy into the gas, which gives us a heating rate 



3.2 x 10 



-28 



c 



10^ 17 s- 



n ergs s 1 cm 3 . 



(46) 



We do not include the effects of heating by soft X-rays, as this is of little importance 
compared to photoelectric heating in low density gas, and is negligible in high density, 
optically thick gas. Other potential heat sources, such as the photoionization of carbon or 
silicon by the UV background, are also insignificant and can be neglected. 

A full list of all of the thermal processes included in our model can be found in Table EJ 



3. Tests 



Tests of the hydrodynami cal and MHD a l gorith ms used in ZEUS are presented in 
Stone fc Normanl (1992a,b) and lHawley fc Stond (119951 ) and therefore are not discussed here. 
Some potential problems relating to the tr eatment of r arefaction waves and adiabatic MHD 
shocks in ZEUS have been pointed out by iFalld (120021 ). but we do not anticipate that these 
problems will invalidate the results presented in this paper. In the case of the rarefaction 
errors, our confidence that they will not significantly affect our results rests on the fact that 
most of the relevant chemistry for H 2 formation occurs in regions of compression while re- 
gions of rarefaction are relatively inactive, and so errors in the treatment of the latter will 
have ve ry little effe ct. As for the shock errors, these vanish if an isothermal equation of state 
is used ( lFalldl2002l ). and since the effective equation of state in our simulations is much closer 
to the isothermal case than the adiabatic case (see § 15.41 below), it is reasonable to expect 
that any errors will be small. 

However, it is necessary to ensure that our modifications to the basic ZEUS code operate 
as intended. To verify the modified code, we have performed tests of the advection of the 
species mass densities in non-reacting flows (discussed in § 13.11 below) and of the solution of 
the chemistry and cooling substep in static gas (discussed in § 13. 2p . 



3.1. Advection tests 



In non-reacting flow, Equation [7] reduces to 

Dp, 



Dt -AV... (47) 
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and so the species mass densities should be advected in precisely the same manner as the 
physical density of the gas. To verify that this is indeed the case, we performed a series of 
tests of the advection alg orithm. Our test suite was based on the advection tests used by 



Stone fc Norman! (Il992al ) to validate the advection algorithms used in ZEUS-2D, a prede- 
cessor of the current ZEUS-MP code. These tests include the advection o f a square pulse 
of high density material by a uniform flow, the classi c Sod shock tube test (ISodlll978l ). and 
various tests taken from IColella fc Woodward! (119841 ) that involve the interaction of strong 
shocks. To use these test problems to verify that our species mass densities are advected 
correctly, we modified them to include two extra field variables, p\ and p 2 , representing ar- 
bitrary chemical species, both of which evolve according to Equation [71 Since we wished to 
simulate non-reacting flow, the chemical creation and destruction terms for these variables 
(Ci, C 2 , Di and D 2 ) were set to zero. The fields were initialized such that pi = p and 
p 2 = 10~ 8 p, the tests were run, and then the resulting values of pi and p 2 were examined. 
Correct operation of the code implies that pi/p and p 2 /p should be conserved. In our tests, 
we found that these ratios were generally conserved to within 0.01%. 



3.2. Chemistry and cooling tests 

Our initial tests of the chemistry and cooling substeps focused on verifying that the 
correct reaction rates and heating and cooling rates were computed by the code, given 
various different sets of input parameters. To do this, we constructed a simple test harness, 
using a mixture of Fortran and Perl, which allowed us to run the main chemistry and cooling 
subroutines separately from the remainder of the hydrodynamical code. We also inserted 
debugging statements into these subroutines which caused them to write out the values 
of the various rates into a log file. The values produced were then compared with those 
produced for the same set of input parameters by an independent test implementation of the 
chemical and cooling rates. Disagreements between the two sets of rates could then be easily 
identified and investigated, and the whole process could be repeated for many different sets 
of input parameters. By eliminating the overhead of running the full hydrodynamical code, 
we were able to quickly explore a wide range of input parameters and verify that the code 
operated correctly in each case. This approach proved especially useful for catching simple 
bugs during the code development process, particularly those caused by typographical errors 
in the software, as these were generally present in either the test implementation or the real 
implementation, but were very unlikely to be present in both. 

We also tested the operation of the chemistry and cooling code at a much higher level 
by performing several sets of simulations of static gas using the full hydrodynamical code. In 
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the first set of simulations, we artificially disabled the operation of the cooling subroutines 
and set A = 0, so that the gas would remain at its initial temperature. We then ran a large 
number of simulations of static gas with different initial temperatures, densities, UV fields 
and dust-to-gas ratios. In each case, we ran the simulation for 3 x 10 15 s, or almost 100 Myr, 
which even at the lowest densities that we examined greatly exceeded the time required for 
the chemical abundances to reach equilibrium. We then compared the value of the H2 and 
H + abundances produced by the simulation with the predicted equilibrium values, which are 
easy to compute for a fixed temperature and density. In every case that we examined, we 
found good agreement between the predicted and the actual values. 

In our second set of simulations, we performed a similar series of tests of the cooling 
routines: we artificially disabled the chemistry, equivalent to setting C» = Dj = in Equa- 
tion [7| above, and solved for the equilibrium temperature of the gas for various different 
initial densities, temperatures, UV fields and H 2 fractions. Again, the resulting values could 
be easily compared with the predicted equilibrium values. Good agreement was again found 
in every case. 

Finally, we performed a set of simulations without artificial constraints on either the 
chemical or the cooling terms. The aim of this set of simulations was t o reproduce the 
phase diagram for interstellar atomic gas calculated by lWolfire et all (120031 ). To do this, we 
ran simulations for a wide range of initial densities, but used the same initial temperature, 
T; = 1000 K, in each case. The other free para meters of our s i mulat ions were set up to 
correspond as closely as possible to those used by IWolfire et al\ ( 120031 ) in computing their 
R = 8.5 kpc Galactrocentric radius model. Specifically, we adopted gas phase carbon and 
oxygen abundances of xq = 1.4 x 10~ 4 and xo = 3.2 x 10~ 4 respectively, a cosmic ray 
ionization rate ( = 1.8 x 10~_ s" 1 and an ultraviolet field strength x — 1-7 (or in other 
words, the standard iDraind (119781 ) field). 



The ph ase diagram we obta ined is plotted in Figure [H We also plot the corresponding 
values from lWolfire et al\ (120031) . taken from their Figure 10. Although our results do not 
match those of lWolfire et al\ precisely, the differences are small and almost c ertainly result 



from the fact that we do not include as wide a range of physical processes as IWolfire et al 



for instance, they model the chemistry of the gas in far more detail than we are able to, 
and so their values for the free electron abundance, and hence for the photoelectric heating 
rate, will be more accurate than ours. The largest difference occurs at n ~ 1 cm -3 , near the 
middle of th e thermally un stable regime, and even here our computed temperatures differ 
from those of lWolfire et all by no more than a factor of two. Given the approximations made 
in our treatment of the chemistry of the ISM, we consider this degree of agreement to be 
acceptable. 
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Fig. 1. — Gas temperature of the ISM as a function of number density, computed for a 
Galactrocentric radius of 8.5 kpc using our m odified version of ZEUS-MP (solid line). For 
comparison, the results of IWolfire et al\ (120031 ) are plotted (dashed line). Note that since our 
modelling of the gas-phase chemistry and gas-grain interactions is considerably less detailed 
than that of Wolfire et al., we do not expect complete agreement. 
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4. Initial conditions 

4.1. Box size and initial number density 

Since our aim in this paper is to model the transition from atomic to molecular gas, 
we have chosen to consider relatively small volumes of gas, which we visualize as being 
small sub-regions within lar g er, gr avit ationally collapsing regi o ns, su ch as those found in 



the simulations of iKravtsovl (120031 ) or ILL Mac Low fc Klessenl (120051 ). For simplicity, and 



also to allow us to use an FFT-based gravity solver, we take our simulation volumes to be 
cubical, with periodic boundary conditions on all sides. For most of our simulations, we used 
boxes of size 20 pc or 40 pc, although in a small number of cases, we used other sizes. The 
particular size used for each individual simulation is summarized in Table HI This table also 
lists the values adopted for the other main adjustable parameters, discussed in more detail 
below. Our choice of box size was guided by two major considerations. The first of these 
was the fact that as we reduce the size of the box, we increase the stability of the gas with 
respect to gravitational collapse. Indeed, if we choose to make the box too small, it will 
contain less than a Jeans mass of gas and so the gas will simply not collapse. The larger that 
we make the box, the more gas it will contain, and the more representative it will be of a 
real gravitationally unstable region. On the other hand, by reducing the size of the box, we 
increase our physical resolution at any given numerical resolution. This is important since 
the physical resolution of our simulations determines the extent to which we can follow the 
gravitational collapse of the gas (see § 14.3.11 below). Obviously, these two considerations 
conflict, but in choosing our box sizes we have done our best to find an appropriate balance 
between them. 

Within the box, we assume that the gas is initially distributed uniformly, with a number 
density ri[. In the majority of our simulations, we take n- x = 100 cm" 3 , as the inferred mea n 



densities of many observed molecular clouds lie near to this value (Mac Low fc Klessenll2004j ) . 



-3 



However, we also explore the effects of reducing rii, examining values as small as 10 cm" 
fsee § I5.5.4D. At a densi ty of 100 cm" 3 , the atomic gas lies well within the CNM regime 



( iWolfire et a/.lll995l . l2003l ). and so the initial temperature of the gas has essentially no impact 
on the results of the simulations, as the gas rapidly cools, reaching thermal equilibrium within 
0.05 Myr. In most of our simulations, we therefore use a rather arbitrary initial temperature 
T; = 1000 K, but as we demonstrate in § 15.41 simulations with T\ = 100 K produce essentially 
identical results for times t <; 0.05 Myr. 

Finally, we note that in all of the simulations described in § El we start with the gas 
at rest. Therefore, in order to trigger gravitational collapse in runs performed using the 
local shielding approximation, it is necessary to introduce small density perturbations to 
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disrupt the symmetry of the initial conditions (as this would otherwise artificially prevent 
gravitational collapse from occurring). These perturbations are not required in runs per- 
formed using our six-ray approximation, as the pressure-driven flows in these runs serve the 
same purpose (see the discussion in § 15.21 below), but we include them nevertheless. To 
produce these perturbations, we select for each grid zone a random number N in the interval 
[—0.5, 0.5], and then perturb the density of that zone by a factor of (1 + 25N), where 5 is an 
adjustable parameter. In most of our simulations, we set 5 = 0.05, although in § 15.5.11 we 
explore the effects of using larger values of S. 



4.2. Initial magnetic field 



Since there is considerable observational evidence for the presence of dynamically sig- 
nificant ma gneti c fields in interstel l ar gas , much of which is discussed in recent reviews by 
Beckl ( 120011 ) and lHeiles fc Crutcherl (120051 ). it is clearly important to take the effects of mag- 
netic fields into account in our simulations. For this reason, most of our simulations involve 
magnetized gas. 

For simplicity, in all of our MHD simulations we assumed that the field was initially 
uniform and oriented parallel to the z-axis of the simulation. The strength of the field was 
a free parameter, and the values used in our various simulations are summarized in Table HI 
Observational determinations of the local magnetic field strength give a typical value of 
6 ± 2/iG, and so we ensured that our fiducial value for the initial magnetic field strength, 
Bifid = 5.85/iG was consistent with this value. 

The reason for our slightly odd choice of fiducial value is that we wanted to ensure that 
in each of our simul ations the mass-to-flux rati o of the gas would be some simple multiple 
of the critical value (INakano fc Nakamural Il978l ) 

1 



ff) 



crit 



(48) 



at which magnetic pressure balances gravity i n an isothermal slab. For a gas cloud in which 
all of the hydrogen is already fully molecular, ICrutcher et all (120041 ) show that the mass-to- 
flux ratio of the cloud can be written in units of this critical value as 



A 



(M/$) 



-21 -^H, 



76xl0~ 21 — ^ (49) 

where Nr 2 is the H2 column density of the cloud in units of cm -2 and B is the strength of 
the magnetic field, in units of /xG. For a fully atomic cloud, the equivalent expression is 



A 



3.8 x 10~ 21 — |i , 
B 



(50) 
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where Nn is column density of atomic hydrogen. For a simulation with an initial atomic 
hydrogen number density of 100 cm -3 and a box size of 40 pc, we have — 1.23 x 10 22 cm -2 , 
and hence A = 46.82/5;, where B x is the initial magnetic field strength. Therefore, in our 
fiducial example in which = B^ = 5.85 /iG, we hav e = 8 . Obse rvati ons of magnetic 



field s trengths in molecular cloud cores, summarized in ICrutcherl (119991 ) and ICrutcher et al. 



(120041 ). find a smaller mean value, A ~ 2, although there is significant scatter around this 
value. There are also many sight-lines for which there are only upper limits on B and hence 
lower limits on A. 

In view of the uncertainty in the appropriate value of A, we ran a number of simulations 
using larger values of the initial field strength, as described in § 15.5.31 i n order to explore the 
extent to which the H 2 formation rate and morphology depend on the field strength. For 
completeness, we also ran a few simulations in which we set B = 0. Although unrealistic in 
the sense that we expect all interstellar gas to be magnetized to some extent, setting B = 
does serve as a guide to the behaviour of the gas in the more realistic case that B is non-zero 
but is too small to significantly affect the gas dynamics on the scales of interest to us. 



4.3. Resolution issues 



4-3.1. The Truelove criterion 



As discussed in § 14.11 one of the important considerations influencing our choice of box 
size was the desire to be able to acc urately follow the gra vitational collapse of gas over as 
large a range of densities as possible. iTruelove et all ( 119971 ) showed that in order to properly 
resolve collapse and avoid artificial fragmentation, it is necessary to resolve the local Jeans 
length of the gas by at least four grid zones. In other words, collapse is resolved only while 



Ax< ^j(p,T), 

where Ax is the width of a single grid zone, and Lj is the Jeans length, given by 



7T 



V2, 



/Tip 



(51) 



(52) 



where c s is the adiabatic sound speed. In Figure El we plot Lj as a function of density 
for gas which is in thermal and chemical equilibrium, with the effects of self- shielding and 
dust shielding calculated using our local shielding approximation and assuming a zone size 
of 0.1 pc. Although somewhat simplified, this model gives a fair representation of how we 
expect the Jeans length to evolve with density in a real simulation. 
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Fig. 2. — Jeans length, Lj, plotted as a function of density, for gas which is in thermal and 
chemical equilibrium. This plot was produced using the same treatment of gas chemistry, 
heating and cooling as was used in our simulations, and using our local approximation to 
treat H 2 self-shielding and dust shielding. 
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It is clear from the figure that it becomes increasingly difficult to resolve Lj as we move 
to higher densities, and so the smaller that we can make Ax, the further we will be able to 
follow the collapse. Now, for a uniform cubical grid, Ax = L/N 1 ^ 3 , where L is the box size 
and N is the total number of grid zones, so we can make Ax small either by increasing iV or 
by decreasing L. Since there are practical limits on how large we can make N, dictated by 
the computational resources required to run the simulation, the easiest way to increase our 
resolution is by decreasing L. However, as noted previously, if we make L too small, then 
the total mass of gas in the box will fall far below the initial Jeans mass, and so the gas will 
be gravitationally stable and will not collapse. 

For any given combination of N and L, we can use Figure [2] to determine the maximum 
number density, n max , at which the Truelove criterion is still satisfied. We have computed 
n max for a range of different numerical resolutions between 64 3 and 512 3 for boxes of size 
20, 40 and 60 pc, and summarize the results in Table [3j The values of n max we obtain vary 
from 490 cm -3 for a 64 3 simulation in a 60 pc box to 5.6 x 10 4 cm -3 for a 512 3 simulation 
in a 20 pc box. Note, however, that these figures are only exact for gas which is in chemical 
equilibrium; if the H2 fraction is out of equilibrium, then this will affect the value of c s , 
through its dependence on 7 and on the mean molecular mass, and so these numbers will 
change slightly. 



4-3.2. The cooling length and the Field length 

Aside from Lj, there are two other important length scales which we would like to 
resolve: the cooling length, L coo i, and the Field length, Lp. A convenient definition of the 
cooling length is the product of the cooling time with the sound speed of the gas: 

-^cool = Cs^cool- (53) 

To properly resolve the thermal state of the flow, particularly features involving strong 
thermal gradients (e.g. shock fronts), we must be able to resolve L coo i- In Figure [31 we plot 
L coo i as a function of density, calculated assuming temperatures which are 1 K (dotted line), 
10 K (dashed line) and 100 K (solid line) greater than the thermal equilibrium temperatureo 

It is clear from the figure that the resolution required to resolve the cooling length of the 
flow is much greater than that required to resolve the Jeans length. This is not unexpected, 
since Lj ~ c B t$, and t coo i <C tg for the conditions of interest to us. However, it does mean 



2 Note that for gas in thermal equilibrium, the cooling rate is zero, by definition, and so t coo \ and L coo \ 
are both formally infinite. 
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Fig. 3. — Cooling length, L coo i, plotted as a function of density, for gas which is 1 K (dot- 
ted line), 10 K (dashed line) or 100 K (solid line) warmer than the thermal equilibrium 
temperature. Note that as L coo \ oc A -1 , the cooling length becomes infinite for gas with a 
temperature which is precisely equal to the equilibrium value. 
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that we cannot reasonably expect to use ZEUS-MP to resolve the cooling length in gas with 
a number density n > 100 cm -3 , unless the temperature of the gas is extremely close to its 
equilibrium value, given the range of values that we use for the box size L. Nevertheless, 
we argue that this does not represent a major problem for the simulations presented in this 
paper, provided that we are careful about the inferences that we draw from them. 

The main consequence of failing to resolve L coo \ is that we overproduce warm gas, 
particularly in post-shock regions, since the width of any such region is clearly constrained 
to be at least one grid zone wide. Therefore, any conclusions that we draw about the 
temperature distribution of the gas should be treated with caution, as we know at the 
outset that we will have too much warm gas, and therefore too little cool, dense gas. The 
effects on the dynamics of the flow are difficult to assess, but it is reasonable to suppose 
that gravitational collapse is inhibited to some degree. Similarly, it is likely that less H 2 is 
produced than would be produced in a higher resolution simulation, since the H2 formation 
rate decreases rapidly with temperature for T > 170 K, and also since the artificially warm 
gas is less dense than it would be if it were able to cool. The net effect of this is that our 
simulations produce less dense gas and less H2 within a given timeframe than would be the 
case if we could resolve L coo i- However, the mass of warm, shocked gas in our simulations 
is never large and so we do not expect our failure to resolve L coo \ to significantly affect our 
conclusions regarding the timescale of H2 formation. 

The other length scale of interest is the Field length, which is the length scale on which 
thermal conductio n stabilizes the growth of thermally unstable density perturbations, and 



which is given by (lFieldlll965l ) 



Mxf- (54) 

where k is the coefficient of thermal conductivity. In Figure HJ we plot the value of Lp 
as a function of density for temperatures that are 1, 10 and 100 K greater than the ther- 
mal equilibrium temperature. To compute these values, we used a value of kh = 2.5 x 
IQ 3 J* 1 / 2 ergs cm " 3 K" 1 s " 1 for the coefficient of thermal conductivity of neutral atomic hydro- 



gen, taken from iParkerl (119531 ). and assumed that k would be of a similar order of magnitude 
in fully molecular hydrogen, and in a mixed gas of atoms and molecules. At the temperatures 
and densities of interest here, thermal conduction by electrons is unimportant and can be 
neglected. 



In a recent paper, iKoyama &: Inutsukal (120041 ) argue that in order to properly resolve 



the dynamics of gas in a thermally bistable medium, it is necessary to resolve Lp with at 
least three grid zones in order to obtain converged results. Now, it is clear from Figure H] that 
it is even harder to resolve Lp than it is to resolve L coo i and so even if we were to include the 
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effects of thermal conduction in our simulations, we could not realistically hope to resolve its 
effects. However, we argue that this does not represent a major problem given the particular 
scenario that we are investigating. The reason is that while gas at the densities studied in our 
simulations may be thermally unstable to isobaric perturbations if its temperature is greater 
than the equilibrium temperature, it is not thermally bistable: it will all cool rapidly to the 
cold phase, rather than forming a two phase medium. Moreover, although isobaric thermal 
instability may operate during this initial cooling phase, amplifying density enhancements 
on scales Lp < L < L crit , where L crit = c s t coo \ is the maximum length scale for an isobaric 
perturbation, the final density will be at most a factor (Ti/T eq ) larger than the initial density, 
where 7] is the initial temperature and T eq the equilibrium temperature of the gas, and the 
resulting density perturbations will not be gravitationally bound. 

To see why, consider an isobaric density perturbation with an initial size L x < L CT[t . 
Before its growth, this perturbation will be unstable to gravitational collapse only if L; > Lj, 
which will be possible only if L crrt > Lj, or in other words if t coo i > tg. However, we know 
that in fact t coo i -C tg for the temperatures and densities of interest to us, and so the initial 
perturbation must be gravitationally stable. Now consider the situation once the gas in the 
perturbation has cooled from 2] to T eq : its density will have increased by a factor of (Ti/T eq ), 
while its local Jeans length will have decreased by the same factor, since Lj oc 
the same time, its linear size will only have decreased by a factor of (Tj/Teq) 1 / 3 , assuming an 
approximately spherical perturbation. Clearly, therefore, the perturbation will be less stable 
against gravitational collapse than it was initially. However, for it to become unstable, its 
initial size must satisfy the following inequality 



T 



2/3 



u > ( -f ) l 3 , m 



and since L ; < L crit , this requires that 



T 



2/3 



L cri t > ( ) (56) 



which in turn will only be satisfied if 



^cool „ / -L eq 



2/3 



*>w ■ (67) 

Now, for T; = 1000 K and T eq ~ 65 K - the values appropriate for fiducial run MS256, 
described in § 15.11 below - we have (T eq /Ti) 2 / 3 ~ 0.16. However, if the initial number 
density is 100 cm -3 , as it is in most of our simulations, then t coo i/tg ~ 0.007. Therefore, 
isobaric perturbations starting with this temperature and density remain gravitationally 
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stable throughout their evolution. A similar result can be derived for the other combinations 
of initial temperature and density examined in this paper. 

If the overdensities created by the isobaric thermal instability are not gravitationally 
unstable and therefore do not collapse further once they have reached T eq , then what actually 
happens to them? As the surrounding gas cools towards T eq , the overdense regions will 
expand in order to stay in pressure equilibrium with the lower density gas, and so their 
density will fall until ultimately they become indistinguishable from the surrounding gas. 
The net effec t will therefore simply be the injection of some additional kinetic energy into 



the gas (as in lKritsuk &: Normanll2002af ). The magnitude of the kinetic energy input depends 
on the fraction of the gas processed through overdense regions, but it will never be larger 
than the initial thermal energy content of the gas. Since the main aim of the simulations 
presented in this paper is to allow us to determine how quickly H2 forms in the absence of 
turbulence, neglecting this energy input is probably justified. 



5. Results 

Before we began the study of H2 formation in turbulent, self-gravitating gas that is 
described in paper II, we first spent some time examining the much simpler case of gas that 
was initially at rest. The main aim of these runs was to determine how quickly H 2 would 
form in the absence of turbulence (§ 15. ip . the role that gravity plays in driving this process 
(§ 15.21) and the morphology of the resulting H2 (§ 15.31) . primarily so that we could later 
compare these results with the results of the turbulent runs discussed in paper II. These 
runs also act as more comprehensive, albeit less quantifiable, tests of the code than the tests 
discussed in § [21 thanks to the relative simplicity of the dynamics, unphysical behaviour is 
much easier to spot here than in the turbulent simulations of paper II. We also examined 
the temperature distribution of the gas and how its maximum and minimum temperatures 
evolve with time (§ I5.4p . and determined how sensitive our results are to variations in input 
parameters such as the box size (§ I5.5.ip . the amplitude of the initial density perturbations 
(§ 15.5.21) . the magnetic field strength (§ 15.5.31) and the mean density of the gas in the box 
(§ I5.5.4p . Parameters for the runs discussed here are listed in Table HI 



5.1. The n 2 formation timescale 



To quantify the rate at which H2 forms in our simulations, there are various quantities 
that we might choose to examine. One of the simplest is the volume averaged molecular 
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fraction, (xh 2 }v, which we can calculate simply by summing up the molecular fraction in 
every grid zone and dividing by the total number of zones, i.e. 

(xu 2 )v = ^^2x R2 (i,j,k), (58) 

where xn 2 {i,j, k) represents the molecular fraction in the zone with coordinates k) and 
N is the total number of zones. However, although (xh 2 )v is very easy to calculate, it 
is of limited use to us, as in an inhomogeneous gas, the volume average will tend to be 
dominated by low density regions, while much of the actual gas resides in high density regions. 
Consequently, a more useful quantity to compute is the mass-weighted mean molecular 
fraction, (xh 2 )m, which is given by 

, E l , J , k Pn 2 ^J,k)AV(t,j, k) 
\Xu 2 )m = , 

where pn 2 (i,j, k) is the mass density of H2 in zone (i, j, k), AV(i,j, k) is the volume of zone 
(i, j, k), Mh is the total mass of hydrogen present in the simulation, and where we sum over 
all grid zones. Finally, we might also look at the total mass of H 2 in the simulation, Mh 2 . 
However, since this can be written in terms of (xh 2 )m as 

M H2 = M h (i H2 )m, (60) 

there is no real benefit to be gained from studying Mh 2 rather than (xh 2 )m- 

To begin our investigation, we first selected a set of parameters that we could treat as 
a fiducial example of collapse from static initial conditions. Our aim here is to discuss the 
results of this example in some detail, and then to explore the effects of varying each of the 
main input parameters by focusing on how the results vary compared to this fiducial case. 
Given the large number of simulations that we have run, this strategy seems to us to be 
more efficient than discussing in detail the results of each individual simulation. 

As mentioned previously, the initial density and temperature used for our fiducial runs 
were = 100 cm -3 and 7]^ = 1000 K respectively. A simple calculation shows that 
with these parameters, the initial Jeans length of the gas is Lj = 47.8 pc, which means that 
if we want the gas in our simulation to be gravitationally unstable at t = 0.0, we would 
need to use a simulation volume of size L > 50 pc. However, since the gas in all of our 
simulations very rapidly cools from 7] ^ to a value close to the equilibrium temperature of 
the gas, it is actually sufficient to choose a value for L that is greater than the value of Lj 
in the cool gas. In the simulations which we performed using our local approximation to 
treat photodissociation and photoheating, the equilibrium temperature is initially the same 
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throughout the simulation volume, and for m = 100 cm -3 has a value T eq ~ 65 K. At 
n = 100 cm -3 and T ~ 65 K, we have Lj ~ 12 pc for fully atomic gas, and so in these 
simulations we require L > 12 pc. In practice, we would like to ensure that we have multiple 
Jeans masses of material in our simulation volume, and so for our fiducial runs we chose to 
set L = 40 pc. In the simulations performed using the six-ray approximation, the equilibrium 
temperature is no longer the same everywhere, as the amount of dust shielding, and hence 
the value of the photoelectric heating rate, varies with position within the simulation volume. 
However, even the largest initial value of T eq in these simulations is smaller than the initial 
value of T eq in the runs performed using the local shielding approximation, and so setting 
L = 40 pc again ensures that multiple Jeans masses of gas will be present. The effect of 
varying L is explored in § 15.5.11 

The other main parameter that we had to specify for our fiducial runs was the strength 
of the initial magnetic field. We chose to set -B^m = 5.85 fiG, for the reasons discussed in 
§ 14.21 For this combination of B, n and L, our initial mass-to-flux ratio (in units of the 
critical value) is A = 8. 

In Figure El we plot the evolution with time of (xh 2 )m for our fiducial parameters for 
six different runs. Runs MS64, MS128 and MS256 were performed using the local shielding 
approximation and differ solely in the numerical resolution used in the runs. Runs MS64-RT, 
MS128-RT and MS256-RT were performed using the six-ray approximation and again differ 
from each other only in terms of numerical resolution. 

It is immediately apparent from Figure [5] that the choice of shielding approximation 
makes a significant difference in the outcome of these uniform density simulations. Far 
more H2 is produced at early times in runs which make use of the six-ray approximation 
than in those using the local shielding approximation. In the six-ray runs, the H 2 fraction 
grows steadily with time, with the gas becoming approximately 50% molecular after lOMyr. 
The H-2 formation rate is independent of the numerical resolution of the simulation until 
t ~ 15 Myr, following which the fastest rate is found in the highest resolution simulation 
(although the difference between the three runs is never great). In the runs performed using 
the local shielding approximation, on the other hand, the growth of the H 2 fraction quickly 
stalls, and the value of (xh 2 )m remains small for at least 20 Myr, particularly in the higher 
resolution simulations. At t > 20 Myr, however, (xh 2 )m suddenly begins to increase rapidly 
in these runs, and by the end of the simulations, the values of (xh 2 )m obtained are starting 
to become comparable with those found in the six-ray runs. Finally, it is clear from Figure [5] 
that the results that we obtain when using the local shielding approximation are resolution 
dependent: the higher the resolution of the simulation, the less H2 is produced. This resolu- 
tion dependence is a consequence of the fact that in the local shielding approximation, the 
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self-shielding factor /shield and dust extinction e~ Td ' 1000 depend explicitly on the physical size 
of the grid zone, Ax, with the result that the equilibrium H 2 fraction, i H2 ,eq, also explicitly 
depends on the zone size. 

In Figure [6], we plot the evolution with time of the maximum gas density, p max , found in 
each of these six simulations. We again see that the results depend on the choice of shielding 
approximation: gas in runs performed using the six-ray approximation collapses roughly 5- 
10 Myr earlier than gas in runs performed using the local shielding approximation. Moreover, 
the two sets of runs display a different sensitivity to the numerical resolution used: in the 
six-ray runs, increasing the resolution causes the collapse to occur slightly earlier, while in 
the other set of runs, increasing the resolution delays the collapse. 

One obvious question to ask is why collapse occurs so much earlier in the six-ray runs 
than in the runs performed using the local shielding approximation. The gas temperatures 
in the former runs are smaller than in the latter (see § 15.41 below) , and so the gas is more 
gravitationally unstable, but the temperature difference is not great and it would be sur- 
prising if this were to be responsible. In fact, the true culprit is something quite different. 
When we use the local shielding approximation, the amount of dust shielding, and hence 
the magnitude of the photoelectric heating rate, is purely a function of the local gas density. 
Since the gas starts with a nearly uniform density distribution, this means that the heating 
rate is the same throughout, meaning that the equilibrium gas temperature is also almost 
the same throughout. In our six-ray runs, on the other hand, this is not the case. Gas 
near the edges of the simulation volume is shielded less than gas at the centre and so is 
heated more. This means that at early times the equilibrium temperature of the gas near 
the edges is higher than that of the gas in the centre. Since the density is still almost the 
same throughout, this results in the development of a pressure gradient, acting towards the 
center of the box. Although the resulting flow velocities are small, as Figure [7] makes clear, 
they are large enough to create non-linear overdensities within about 10 Myr, which subse- 
quently merge and collapse on the free-fall timescale of 5 Myr. In comparison, the density 
inhomogeneities present in runs MS64, MS128 and MS256 at t — 10 Myr are very much 
smaller and so gravitational collapse is delayed for a much longer period. 

An interesting consequence of this is that H 2 formation can actually occur faster when an 
ultraviolet field is present then when one is absent, as we show in Figure [HJ Of course, this is 
not a new discovery: what we are seeing is essentially a form of radiation-driven implosion, 
which is a phenomena that has been discussed by a number of previous authors, albeit 
primaril y in the context of driving by photoioniza t ion, rather than by photoelectric heating 
(see e.g. 
2003h . 
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Given the large disparity between the two sets of results, which set of results should 
we believe? In other words, which of the two approximations that we have used does better 
at capturing the true behaviour of the gas? From our discussion in § 12.2.11 we know that 
the main disadvantages of the six-ray approximation are its lack of angular resolution and 
its insensitivity to the effects of velocity differences along the line of sight. However, at 
early times in these simulations, the gas distribution is nearly uniform and the velocities are 
small, and so this approximation should give very accurate results. On the other hand, we 
know that the local approximation underestimates the true amount of shielding, and from 
our results here it is apparent that we underestimate the shielding by quite a large factor, 
enough to significantly alter the outcome of the simulations. We can therefore conclude 
that this approximation does not work well for treating H2 photodissociation in this case. 
Of course, this does not come as a great surprise, as the physical conditions which best 
motivate the use of a local approximation - large variations in the gas velocity from zone to 
zone, and significant density inhomogeneities - are not present at early times in these runs. 
At late times, once a substantial density inhomogeneity has developed, the local shielding 
approximation does a much better job of modelling the photodissociation. Moreover, as 
we will see in paper II, the local shielding approximation works far better for treating H 2 
photodissociation in supersonically turbulent flow than it does here. 

As larger density inhomogeneities develop, and in particular as large flow velocities 
develop, the accuracy of our six-ray approximation degrades. However, Figures [6] and [7] 
demonstrate that the flow remains subsonic and the density inhomogeneities remain small 
at t < 10 Myr, and so the six-ray approximation should remain accurate throughout this 
period. Moreover, at later times enough H2 is present to provide effective self-shielding for 
almost all of the gas, and so we would not expect our results at late times to be particularly 
sensitive to the growing inaccuracy of our approximation. 

It is also clear that when we use the six-ray approximation, we do a much better job 
of modelling the dynamics of the flow, as we recover the large scale pressure gradient that 
is missed by the local approximation. It should be noted, however, that one reason that we 
find so much difference between the dynamics of the two sets of simulations is our choice 
of initial conditions. Because the gas is initially at rest, any small velocities, such as those 
produced by the large-scale pressure gradient, are significant and can have a large effect on 
the outcome. In simulations where the gas does not start at rest (such as the turbulent 
models presented in paper II), we would expect to see far less difference in the dynamics. 

In view of the higher accuracy of the results from our six-ray runs, we will focus on 
these in particular in the remainder of this paper. However, for the purposes of compar- 
ison we include some discussion of results from runs performed using the local shielding 
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approximation 

Finally, it should be noted that even in the six-ray runs, it takes roughly lOMyr to form 
a significant amount of H 2 . This is the most important result of the simulations reported on 
in this paper, as this timescale is greater than a gravitational free-fall time and is significantly 
longer than the timescales derived for the turbulent cloud models studied in paper II. 

5.2. The role of gravity 

From the comparison of Figures [5J and El it is clear that the growth in (o;h 2 )m at late 
times in runs MS64, MS128 and MS256 is being driven by the gravitational collapse of the 
gas: the increase in density caused by the collapse increases both the H2 formation rate and 
the amount of shielding, allowing (xh 2 )m to rise rapidly. Similarly, the growth of (jh 2 }m at 
late times in the six-ray runs appears to be accelerated by the collapse of the gas. 

To confirm this, we performed a set of runs in which we did not include the effects 
of self-gravity. These runs, designated MS64-ng, MS128-ng, MS256-ng and MS256-RT-ng, 
used the same input parameters as runs MS64, MS128, MS256 and MS256-RT respectively. 
The evolution of (x H2 >m in runs MS256, MS256-ng, MS256-RT and MS256-RT-ng is shown 
in Figure El At t < 15 Myr in runs MS256-RT and MS256-RT-ng, and at t < 27 Myr in 
runs MS256 and MS256-ng, the results of the simulations agree, demonstrating that at early 
times the effects of self-gravity are unimportant. However, it is also clear that at later times, 
the simulation results diverge, with gravitational collapse providing a distinct boost to the 
H 2 formation rate in the self-gravitating simulations. 

Since gravity plays an important role in the outcome of our simulations, it is also natural 
to investigate how well we are resolving the gravitational collapse of the gas. Specifically, 
we would like to know at what point the Truelove criterion is first violated and how much 
of the H 2 formation that we see in our runs occurs before this point. 

Using the results of our simulations, it is relatively easy to determine when the Truelove 
criterion is first violated for any given model, and what the value of (xh 2 )m is at that time. 
The values we obtain for each of our runs are summarized in Table [51 In runs performed 
using the local shielding approximation, most H 2 formation occurs after the run has violated 
the Truelove criterion. On the other hand, in runs that use the six-ray approximation, the 
majority of the H 2 forms before the Truelove criterion is violated. 

Of course, the fact that we no longer properly resolve gravitational collapse in dense gas 
once the Truelove criterion is violated does not necessarily invalidate all of the subsequent 
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results from our simulations. For instance, if only a few percent of the total H2 were found in 
unresolved regions, then the fact that we are no longer able to resolve these regions properly 
would be likely to have only a minor effect on (xh 2 )m- Conversely, if most of the H 2 were in 
unresolved regions, then our results for (x H2 )m would be trustworthy only up to the point 
at which the Truelove criterion was violated. In order to quantify how much gas and how 
much H 2 ends up in unresolved regions in our simulations, we examined intermediate output 
dumps of density, internal energy and H 2 fraction from each of our fiducial runs, and for 
each dump determined which zones, if any, were unresolved. We then computed f res , the 
fraction of the total gas mass in resolved regions, and / r es,H 2 the fraction of the total H 2 
mass in the same resolved regions, for each output time. The resulting values of / rcs and 
/res,H 2 are plotted in figures [TUh and [TUb respectively. It is clear from these figures that once 
the Truelove criterion is violated, the majority of the gas and of the H 2 soon ends up being 
located in unresolved regions, suggesting that our simulation results are only trustworthy up 
to this point. 

We obtain similar results if we perform the same analysis for any of the other runs in 
Table H] which violate the Truelove criterion. We should therefore not place too much weight 
on results coming from very late times in any of these simulations. Nevertheless, it should 
be clear that these concerns do not affect the basic timescale for H 2 formation that we find 
from our simulations and so do not undermine our main results. 
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Fig. 4. — Thermal conduction length, or Field length, L F , calculated assuming temperatures 
which are IK (dotted line), 10 K (dashed line) and 100 K (solid line) greater than the thermal 
equilibrium temperature. 
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Fig. 5. — Time evolution of (x H2 )m in several sets of runs. All of the simulations were 
performed using our fiducial set of initial conditions. The three lines in the bottom right 
of the plot show the results from runs MS64 (dot-dashed line), MS 128 (dashed line) and 
MS256 (solid line), performed using the local shielding approximation. The three lines in 
the upper left indicate the results of runs MS64-RT (dot-dashed line), MS128-RT (dashed 
line) and MS256-RT (solid line), performed using the six-ray shielding approximation. The 
digits in the run name indicate the numerical resolution; i.e. runs MS64 and MS64-RT were 
both performed with 64 3 zones resolution etc. 
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Fig. 6. — Time evolution of the maximum gas density, p max - The three lines on the left-hand 
side of the plot correspond to runs MS64-RT (dot-dashed line), MS128-RT (dashed line) and 
MS256-RT (solid line); the lines on the right-hand side correspond to runs MS64 (dot-dashed 
line), MS128 (dashed line) and MS256 (solid line). 
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Fig. 7. — Time evolution of the maximum velocity in 256 3 runs MS256-RT (solid line) and 
MS256 (dot-dashed line). The evolution of the sound speed in the densest gas in run MS256 
is indicated by the dotted line; comparable results are found for the sound speed in run 
MS256-RT, but we omit them here for clarity. The larger maximum velocities found in the 
six-ray run at early times are caused by pressure waves driven in from the boundaries by 
photoelectric heating there. 
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Fig. 8. — Time evolution of (xh 2 )m in 256 3 runs MS256-nr (solid line), MS256-RT (dashed 
line) and MS256 (dot-dashed line). Run MS256-nr was performed with the strength of the 
radiation field set to zero. The fact that more H 2 forms in run MS256-RT than in run 
MS256-nr is attributable to the compression of the cloud that occurs in the former case due 
to the temperature and pressure gradients set up by the non-uniform heating of the gas. 
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Fig. 9. — Time evolution of (£h 2 )m m 256 3 runs performed with and without self-gravity. 
The two lines on the left-hand side of the plot indicate the results of runs MS256-RT and 
MS256-RT-ng (solid and dashed lines respectively). The two lines on the right-hand side 
indicate the results of runs MS256 (solid line) and MS256-ng (dashed line). Runs MS256-RT 
and MS256 include the effects of self-gravity, while runs MS256-RT-ng and MS256-ng do 
not. 
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Fig. 10. — (a) Fraction of the total gas mass situated in resolved regions (i.e. in zones which 
satisfy the Truelove criterion) plotted as a function of time for runs MS64-RT (left-hand 
dot-dashed line), MS128-RT (left-hand dashed line), MS256-RT (left-hand solid line), MS64 
(right-hand dot-dashed line), MS128 (right-hand dashed line) and MS256 (right-hand solid 
line), (b) As (a), but for the fraction of the total mass of molecular gas. 
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5.3. H2 morphology 

As noted above, even in our highest resolution fiducial simulations, most of the gas is 
located in unresolved regions at the end of the run. It is therefore not particularly enlight- 
ening to examine the morphology of either the gas or the H 2 at the end of the run, as both 
will be inaccurate in ways not easily quantified. Moreover, in a real molecular cloud, we 
would expect star formation to occur rapidly once gas has collapsed to high densities, and 
since we do not include any feedback in these simulations from effects such as protostellar 
outflows, our simulations are missing some of the physics necessary to accurately model the 
morphology at late times. 

On the other hand, what we can do with some degree of confidence is to examine the 
morphology of the gas at a time shortly before the Truelove criterion is first violated, when 
star formation has presumably yet to occur. The morphology of the cloud at this time 
depends upon our choice of shielding approximation. In the runs performed with the six-ray 
approximation, a thick slab of gas has formed, oriented perpendicularly to the magnetic 
field. The density within most of this slab is only slightly elevated over the initial density 
of the gas, but the thin layer of gas bounding the slab shows a large overdensity. Plots of 
the gas density in the x-z and x-y planes in run MS256-RT at time t = 17.4 Myr are shown 
in Figure [TTJ The gas in this simulation is gravitationally collapsing in a direction parallel 
to the magnetic field lines, with the result that shortly after the time of this output dump, 
the two large overdensities merge as the gas forms a thin dense sheet located at z = (see 
Figure [T2|) . It is reas onable to exp ect that this sheet would then fragment into a number of 



:xp 

m 



filaments and cores (|Larsonl 11985ft . but by this point in our simulation, the sheet itself was 



unresolved and so we were unable to follow the further dynamical evolution of the simulated 
cloud. 

The morphology of the gas in run MS256 at an output time shortly before the Truelove 
criterion is violated is also sheet-like, as we can see from Figure [131 However, the width of 
the sheet in run MS256 is much smaller, and the density cross-section is rather different, 
with the maximum density being found at the midplane, rather than at the top and bottom 
edges. As we have already discussed, these differences stem from the differing dynamics of 
the flow: in this run, photoelectric heating by the ultraviolet background is initially uniform, 
and so there is no large scale temperature or pressure gradient. The fact that in both runs 
gravitational collapse produces a sheet-like structure is a consequence of the presence of the 
magnetic field, which strongly suppresses gas flow perpendicular to the field, but which does 
not affect flow in the z direction, parallel to the field. 



The spatial distribution of H 2 in these simulations is also of interest. In Figure 
we show how the H 2 fraction varies with position in run MS256-RT. Figure [TB] gives the 
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Fig. 11. — (a) Slice in the x-z plane through the density field in 256 3 zone run MS256-RT 
at a time t = 17.4 Myr. The slice is centered on the midpoint of the simulation volume, (b) 
As (a), but for a slice in the x-y plane. Note the change in density scale compared to figure 
(a). 



-50- 



corresponding results for run MS256. In both runs, there is an obvious correlation between 
the gas density and the H 2 fraction, with the largest molecular fractions being found in the 
densest gas. However, this correlation is much stronger in run MS256 than in run MS256- 
RT, because in the former run, the amount of shielding is purely a function of the local 
gas density, so H 2 both forms faster and is photodissociated more slowly in denser gas. In 
run MS256-RT, the correlation between local gas density and photodissociation rate is much 
weaker (although not completely absent). 

To better quantify the relationship between gas density and H2 fraction, we examined 
how the mean H2 fraction varied as a function of n in each simulation. To do this, we 
computed xn 2 and n for each grid zone, and then binned the data by number density, using 
bins of width 0.05 dex. We then computed the mean and standard deviation for x^ 2 for 
each bin. The resulting values are plotted is Figure [T6k for run MS256-RT and Figure [T6b 
for run MS256. Although the mean values that we compute here are volume weighted, we 
would not expect the mass weighted values to differ greatly, since the narrow width of our 
density bins means that there is little variation in the gas mass from zone to zone within a 
given bin. 

An immediately obvious feature of Figure [T6a is the discontinuity at n = 100 cm -3 . 
This feature is a result of the fact that some of the H2 that forms in the overdense regions 
bounding the collapsing slab is left behind by the collapse, and so finds itself ultimately in 
a region with n < 100 cm" 3 . In other words, most of the gas at n < 100 cm -3 did not form 
in situ, but instead formed at higher densities and has been transported to lower densities. 
On the other hand, gas with n > 100 cm -3 resides in the dense slab, and much of this gas 
has yet to be greatly affected by the collapse (as is apparent from Figure [HI). 

We do not see a comparable feature in Figure [TdT d. This is not unexpected, as the dy- 
namics of the flow are quite different in this case, with collapse happening far more slowly. 
Moreover, the local shielding approximation used in run MS256 leads naturally to a tight 
relationship between density and H 2 fraction at low densities, as the gas is close to photodis- 
sociation equilibrium. This can be seen clearly in Figure [T6b if one compare the simulation 
results with the curve indicating how £H 2 ,eq varies with density. The latter curve was con- 
structed by first computing £H 2 ,eq for each grid zone (assuming shielding to be described by 
our local shielding approximation) and then binning and averaging these values using the 
same procedure as for xh 2 . Note that in run MS256-RT we expect £H 2 ,eq to vary with both 
density and with position within the simulation volume, and so we have not constructed a 
similar curve for this run. 

It is also of interest to ask how the H2 mass is distributed as a function of density. The 
plots given above show how the H2 fraction varies with density, but contain no information 
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on how much gas is to be found at any given density. This information is given by the 
mass- weighted density probability distribution function (PDF), which we have plotted for 
runs MS256-RT and MS256 in Figures [T7h and [TTb respectively. 

Figure [T7| demonstrates that in both runs, much of the gas remains close to the initial 
density. In run MS256-RT there is also a substantial amount of gas at n > 1000 cm -3 , 
corresponding to gas associated with the overdensities bounding the collapsing slab. In run 
MS256, this feature is absent, but there is nevertheless a high density tail, extending up to 
n ~ 4000 cm -3 . In run MS256-RT, about 30% of the total H 2 mass is associated with the 
relatively unperturbed gas within the slab, with the other 70% being associated with the 
overdense, collapsing gas (see Figure [T8l) . On the other hand, in run MS256, most of the H2 
is located at densities close to the initial gas density; indeed, only half of the total amount 
of H 2 is found at n > 100 cm -3 . 

However, as the gas continues to gravitationally collapse, we expect both the PDF and 
the cumulative mass distribution of H2 to alter greatly, as we know already that most of the 
mass and most of the H 2 will ultimately be located in dense, unresolved gas. This is borne 
out by the results from the end of the simulations which are also plotted in Figures [T7] and 
M for runs MS256-RT and MS256. 

5.4. Gas temperature: evolution and distribution 

As we have previously noted, the cooling time of the gas in our simulations is much 
shorter than the dynamical timescale. Therefore, the gas very quickly cools to the thermal 
equilibrium temperature (which for a gas density of 100 cm -3 is approximately 65 K, if we use 
our local shielding approximation) regardless of the initial temperature of the cloud. This 
is clearly illustrated in Figure [T9], where we plot the evolution with time of the minimum 
and maximum gas temperatures T min and T max in run MS256 (which has an initial gas 
temperature 7] = 1000 K) and run MS256-T100 (which has T ; = 100K but is otherwise 
identical to run MS256). Following an initial period of cooling that lasts for approximately 
0.05 Myr, the behaviour of the two runs becomes essentially identical. We see a similarly 
rapid initial phase of cooling in runs performed using the six-ray shielding approximation. 

In run MS256, the temperature distribution following the initial cooling phase is almost 
uniform, and remains so until the runaway gravitational collapse of the gas begins at t ~ 
20 Myr. On the other hand, in run MS256-RT we see that even at early times there is a 
temperature difference of approximately 20 K between T m [ n and T max . This temperature 
differential is a result of the fact that gas at the center of the box is more shielded from 
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photoelectric heating than gas at the edge, and, as we have already discussed, this gives rise 
to a pressure gradient that drives the subsequent gas flow. 

In both runs, we find tight correlations between the gas density and the temperature, 
as we demonstrate in Figure [20j In run MS256-RT, we see that there is a discontinuity at 
n = 100 cm -3 . Gas at densities n > 100 cm -3 is located within the collapsing slab, and so 
is shielded from the ultraviolet background by the gas in the overdense regions that bound 
the slab. On the other hand, gas at n < 100cm -3 is located either above or below the slab 
and so only receives the benefit of this shielding in one direction. Since the contribution 
from the lightly shielded direction dominates, the result is that the photoelectric heating 
rate is considerably higher in this gas than in the dense gas in the slab, and so consequently 
T eq is also higher. No comparable effect is seen in run MS256, as this is a consequence 
of the non-local nature of the dust shielding and so is not captured by the local shielding 
approximation. 

We have also examined whether one can usefully describe the behaviour of gas in these 
simulations using a polytropic equation of state. In other words, if we write the gas pressure 
clS cL polynomial function of density: 

p = Kp™, (61) 

then what functional form must 7 e g have if this relation is to accurately describe the gas? If 
7 e ff is a constant, independent of the density, then the gas is a simple polytrope. Even if j e g 
is not constant, however, this description can be useful provided that 7 e fr varies smoothly 
and simply with p. 

To investigate how / y c g varies as a function of density at late times in each run, we 
computed the thermal pressure p for each grid zone and then binned the data with density 
just as we did for the temperature above. We then used the resulting curve, together with 
the fact that 

7cfr = -t] — (62) 
dmp 

to compute 7 e fj as a function of density. The values we obtained for runs MS256-RT and 
MS256 are plotted in Figures l21h and 121b respectively. 

Figure l2Ta demonstrates clearly that the gas in run MS256-RT is not well described by 
a polytropic equation of state. The prominent feature at n ~ 100 cm -3 is a consequence of 
the discontinuity in the temperature curve at this density and should perhaps be disregarded 
(since one might still hope to use a polytropic description for the gas within and without 
the slab, even if the polytropic index for one differs from that for the other). However, from 
the figure it is clear that even at densities n ^> 100 cm -3 and n <C 100 cm -3 a polytropic 
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equation of state is not really appropriate, as there is considerable variation in 7 e g with 
n. Nevertheless, it is clear that the equation of state of the gas is softer than isothermal, 
considerably so at low densities. 

In run MS256, the polytropic approximation fares much better. Although the gas is 
not a polytrope, since 7 e fr varies with density, the rate of change of 7 e fr is typically small. 
For the range of gas densities covered by our simulations, 7 e fj ~ 0.7-0.8, and so for many 
applications, treating the gas as a simple polytrope with a constant polytropic exponent that 
is within this range of values may be a reasonable approximation. We hypothesize that the 
polytropic description does much better in this case because when we use the local shielding 
approximation, all of the heating and cooling terms in the energy equation become smooth 
functions of density. This is not the case in the six-ray runs, as the non-local shielding 
disrupts the simple relationship between density and photoelectric heating rate. 

It is also interesting to examine how the values of 7 e fr derived here compare with previous 
suggestions in the li terature. Our values are significantly smaller than those derived by 



Spaans fc Silkl (120001 ) from their numerical models of solar metallicity gas with this range of 



densities. However, the values found in run M S256 are only slightly smaller than the value 



of 7 e ff = 0.725 derived by lJappsen et al\ (120051 ) from a synthesis of a variety of observational 



and theoretical sources. In run MS256-RT, the value of 7 e g at low densities is considerably 
smaller still, but at higher densities there is better agreement. Nevertheless, it is important 
to bear in mind that our treatment of the gas chemistry remains highly approximate, and 
that this, plus the crude nature of our shielding approximations, will limit the accuracy with 
which we can model the gas temperature. It is the refore unclear h ow se rio usly we should 



take t he differences between our results and those of ISpaans fc Silkl (120001 ) or lJappsen et al. 



(120051 ) . It would be interesting to revisit this point in the future with simulations that 
include far more of the relevant carbon and oxygen chemistry and a better treatment of the 
photoelectric heating. 

Finally, we have investigated how the current uncertainty in the value of the cosmic 
ray ionization rate, discussed previously in § 12.21 affects the temperature evolution of the 
gas in our simulations. To do this, we performed a run, designated MS256-CR, with the 
same input parameters as run MS256, but taking a large value of £ = 10~ 15 s^ 1 for the 
cos mic ray ionization r ate. This value lies at the high end of current determinations (see 



e.g. iMcCall et a/.ll2003l ) and so this run and run MS256 should bracket the true behaviour. 
We compare the evolution of T min and T max in these two runs in Figure [221 We see from 
the figure that prior to the onset of runaway gravitational collapse, the effect of the higher 
cosmic ray ionization rate is to increase both T min and T max by about 20 K. However, once 
collapse begins, the difference between the runs is quickly erased. 
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The higher temperatures resulting from the higher ionization rate have some influence 
on the H 2 formation rate. At early times, the net effect is small: the higher gas temperature 
leads to a slightly higher rate coefficient for H 2 formation, and so to a higher H 2 abundance, 
but the difference is at the level of about 10%. At late times, once gravitational collapse is 
well underway, the numbers reported in Table [5] suggest a rather more substantial difference 
between the runs. However, most of this difference is illusory, and is a consequence of the 
fact that the higher temperature of the gas in run MS256-CR is matched by a higher Jeans 
mass, allowing us to resolve the collapse for about 0.5 Myr longer. We therefore resolve more 
of the H-2 formation that occurs during this run. 

We also briefly investigated the effects of using a higher value of ( in combination with 
our six-ray approximation. However, low resolution test runs suggested that the differences 
were no larger than the differences between runs MS256 and MS256-CR, so we did not pursue 
this line of investigation further. 

5.5. Sensitivity to variations of the input parameters 

Having explored in some detail the formation and distribution of H 2 in our standard 
runs, and the thermal behaviour of the gas, it is now time to turn our attention to examining 
what happens when our input parameters are varied from their standard values. Therefore, in 
§ 15. 5. 1145. 5.41 below, we examine the effects of varying the box size L, the initial perturbation 
amplitude 5, the initial magnetic field strength B\ and the initial density n\, while holding 
the other parameters constant. 

5.5.1. Box size 

By varying the size of the box while keeping the density of gas within it fixed, we 
can alter the number of Jeans masses of gas which the box contains and so make the gas 
either more stable against gravitational collapse (if we decrease L) or less stable against 
collapse (if we increase L). If the results of our fiducial simulations are to be considered 
to be properly representative of the behaviour of real gas, then we should ensure that they 
are insensitive to small changes in L. We therefore performed several additional simulations 
with our standard set of input parameters, but using different values for L. Specifically, we 
performed four additional 256 3 runs using the local shielding approximation, with L = 20, 
30, 50 and 60 pc respectively, and two additional runs using the six-ray approximation, with 
L = 20 and 60 pc respectively. The results of these runs are plotted in Figure [231 along with 
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the results from runs MS256 and MS256-RT for comparison. 

We see that in the runs using the local shielding approximation, increasing L at fixed 
numerical resolution increases the value of (xh 2 )m at all times. This is easy to understand 
as in these runs, the amount of shielding is proportional to the grid zone size Ax, and if 
we increase L while keeping the numerical resolution fixed, then we increase Ax. It is also 
apparent that in these runs the rapid increase in (xh 2 )m that is caused by the runaway 
collapse of the gas occurs earlier in runs with a larger value of L. Figure I2"3"b demonstrates 
that this is because runaway gravitational collapse occurs at progressively earlier times as 
L is increased and at later times as L is decreased. This is a straightforward consequence 
of the presence of an increased number of Jeans masses in the larger simulations compared 
to the smaller ones. However, it should be noted that the changes to L change the collapse 
time of the gas by no more than 20%. Therefore, while magnetic and thermal pressure are 
clearly playing some role in retarding the collapse of the gas in our fiducial simulations, the 
effect is not so large as to render the results of these runs unrepresentative of the general 
case, particularly given the level of approximation to which we are working. 

In the runs using the six-ray approximation, we see much less sensitivity to the box 
size, which again is not unexpected, given that the shielding in these runs is not directly 
dependent on Ax. The main effect that is apparent is that H2 forms slightly more efficiently 
at early times and considerably more efficiently at late times in runs with a smaller L. This 
appears to be a consequence of the pressure-driven dynamics of the flow, which we have 
already discussed in some detail in previous sections. In runs with a smaller L it takes less 
time for the overdense gas layers to propagate to the center of the box, and so therefore it 
also takes less time for the gas to reach a state of runaway gravitational collapse, as can 
be seen clearly in Figure [23b . The more rapid density evolution in runs with smaller L 
leads to a more rapid growth of H 2 . Nevertheless, the difference between the runs remains 
relatively small, suggesting again that our fiducial simulation is adequately representative of 
the general case. 

5.5.2. Initial perturbation amplitude 

Since we use periodic boundary conditions in our simulations, a perfectly uniform dis- 
tribution of gas will be in dynamical equilibrium (albeit an unstable equilibrium on scales 
larger than the Jeans length) and so will not collapse. To provoke collapse, it is neces- 
sary to perturb the distribution. In our simulations of turbulent gas, described in paper 
II, large perturbations are rapidly created by the motion of the gas. Similarly, in the runs 
described here which used the six-ray shielding approximation, the pressure gradient caused 
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by the non-uniform photoelectric heating drives a flow that creates density perturbations 
large enough to trigger gravitational collapse. However, in the runs we performed that used 
the local shielding approximation, it was necessary to add some form of initial perturbation 
by hand. Our technique for doing so has already been described, but it is important to 
understand to what extent our results depend on our choice of 5, the maximum amplitude of 
the initial random perturbations. Specifically, we would like to know whether the long delay 
before runaway gravitational collapse begins in earnest in our simulations is a consequence 
of our choice of a small value for 5. Therefore, we have run several simulations in which we 
have varied S while keeping the other parameters fixed at their fiducial values. 

The growth of (xh 2 )m with time in 256 3 zone runs with 5 = 0.1, 0.5 and 1.00, as well 
as in run MS256 is shown in Figure [24] Clearly, increasing 5 does decrease the time that 
elapses before runaway collapse occurs, as can also be seen from the results summarized 
in Table [5j Nevertheless, the time required to form significant quantities of H2 remains 
long even in the 5 = 1.00 case, and while we could reduce the timescale even further by 
considering even larger inhomogeneities, at this point we would essentially be considering 
a collection of smaller, denser clouds, rather than a single cloud with a well-defined initial 
density. Therefore, although the timescale for collapse in these runs does depend on the 
initial density structure of our cloud, we can nevertheless be certain that it is of the order 
of 20 Myr or more for an approximately uniform cloud initially at rest. 

We also investigated the sensitivity of runs that use our six-ray approximation to varia- 
tions in 5, since these runs are also seeded with small perturbations to the initially uniform 
density. As expected, these runs show essentially no sensitivity to 5: the aforementioned 
pressure- driven effects dominate. 

5. 5. 3. Initial magnetic field strength 

In order to explore how our results depend on the strength of the initial magnetic field, 
we performed several sets of runs in which we varied the strength of the field but kept all of 
the other input parameters the same as in our fiducial runs. Specifically, we performed two 
additional runs using the six-ray shielding approximation: one with an initial field strength 
-Bi,M — 23.4 /iG, corresponding to a mass-to-flux ratio of two, in units of the critical value, 
and one with a field strength of zero, i.e. a purely hydrodynamical run. The evolution with 
time of (xh 2 )m in these runs is plotted in Figure [25k. along with the corresponding values 
from run MS256-RT. We see that at t > 14 Myr, the H2 fraction rises more rapidly in our 
hydrodynamical run than in our MHD runs, and also that the behaviour of the latter runs 
is essentially indistinguishable. Figure [23b . which shows the evolution with time of p max i n 
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these three runs, demonstrates that difference in behaviour between the run with B = and 
the runs with B ^ is caused by the fact that gravitational collapse occurs more rapidly 
in the former run than in the latter runs. This is a simple consequence of the fact that in 
the absence of the magnetic field, gas can flow freely in every direction, allowing the cloud 
to collapse along all three of its axes at once, while when a magnetic field is present, the gas 
flow is channeled primarily along the field lines, leading to collapse along only one axis. This 
also explains why our results are not particularly sensitive to the strength of the field, when 
one is present, as this will have no effect on the velocity of the flow parallel to the field lines. 

We also examined the sensitivity of the results from the local shielding approximation 
simulations to the value of B\, using values of 0.0, 11.7 /zG and 23.4 /zG, and found very 
similar results: collapse occurs faster and H 2 forms more quickly in the complete absence of 
a magnetic field, but if a field is present then the timescales for either are not very sensitive 
to its strength. 

5.5.4- Initial density 

A final interesting topic to examine is the sensitivity of our results to the initial density 
of the gas. For instance, if we reduce the initial number density rii from its fiducial value of 
100 cm -3 , which is on the high side for the CNM, to 30 cm -3 or 10 cm -3 , what effect will this 
have on the timescale for H 2 formation? Obviously, if we reduced n ; but kept L constant, 
then we would reduce the amount of gas in our simulation volume and thereby alter how well 
the gas could resist gravitational collapse, just as we did when we reduced L while keeping 
n\ constant. Therefore, to examine the effect of changing the density without significantly 
affecting the gravitational stability of the gas, we must increase L at the same time that we 
decrease n\. Similarly, if we decrease n\ and increase L, we must also alter B\ if we wish 
to keep the mass-to-flux ratio constant, and hence ensure that the field cannot completely 
prevent the gas from collapsing. 

We therefore performed three simulations with lower initial densities than in our fiducial 
simulations. In run MS256-RT-nlO, we set n\ = 10 cm -3 and used our six-ray shielding 
approximation. In runs MS256-nlO and MS256-n30, we set n\ = 10 and 30 cm -3 respectively, 
and used our local shielding approximation. In all three cases, the box size and magnetic field 
strength were adjusted so as to keep the number of Jeans masses in the simulation volume 
and the mass-to-flux ratio approximately the same as in our fiducial runs (see Table H] for 
the values used). The evolution of (xh 2 )m in these simulations, as well as in runs MS256 and 
MS256-RT for comparison, is plotted in Figure [26k . In Figure [26b . we show a similar plot 
for p max 
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Comparison of the evolution of (o;h 2 )m in the various runs (see Figure [26k) demonstrates 
that the reduction in n; has a huge effect on the evolution of (xh 2 )m at early times in runs 
performed using the local shielding approximation. This is a consequence of the dramatic 
reduction in £H 2 ,eq in these runs, which decreases from ~ 0.1 in run MS256 to ~ 4 x 10~ 5 
in run MS256-n30 and to ~ 7 x 10~ 6 in run MS256-nlO. Substantial quantities of H2 are 
produced in run MS256-n30 only once the gravitational collapse of the gas is well underway 
and the density in the collapsing gas has become high enough to allow for effective self- 
shielding on small scales. Since the reduction in n\ from 100 cm -3 to 30 cm" 3 more than 
doubles the time required for the gas to collapse, this means that we do not see significant 
H2 formation before t ~ 40 Myr. In run MS256-nlO, substantial quantities of H2 are never 
produced, even after 100 Myr (although we suspect that if we were to run the simulation for 
considerably longer, we would see significant H 2 formation, as the gas in this run has only 
just begun to collapse at the point when the simulation is ended). 

In run MS256-RT-nlO, the effect of the reduced density is less pronounced, but is signifi- 
cant nonetheless. The reduction of n\ by an order of magnitude leads to a reduction in (xh 2 )m 
by roughly a factor of seven at early times. The reduction in n x again leads to an increase 
in the time required for gravitational collapse to occur, and the gas in run MT256-RT-nlO 
does not become dominated by H2 until collapse occurs at t ~ 50 Myr. 

In practice, of course, we do not expect the gravitational collapse of gas clouds with 
mean densities of 10 or 30 cm -3 to take quite as long as these simulations suggest, as the 
use of periodic boundary conditions to treat real clouds with these densities, and with the 
sizes considered here, is unlikely to be a good approximation - there is simply too much 
inhomogeneity in the ISM on these scalesj. Nevertheless, even if we were to disable the 
periodic boundary conditions and to treat the g isolated cloud in free-fall collapse 

(i.e. if we could ignore its thermal and magnetic energy content), then we would still expect 
gravitational collapse, together with the associated H 2 formation, to occur on a timescale of 
the order of 10 Myr or more. We show in paper II that the results of turbulent simulations 
with the same mean density are dramatically different. 

6. Summary 

In this paper, we have discussed in detail how we have modified the ZEUS-MP hydro- 
dynamical code to allow us to simulate the formation and destruction of H2 in the ISM. 



3 For instance, the scale height of molecular gas in the disk of the Mil ky Way at the solar Galactrocentric 
radius is estimated to be only 65 pc ( Sanders. Solomon fc Scovillel 1984 ). 
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We have used two different approaches for modeling the shielding of H2 against photodis- 
sociation. The first of these is a local shielding approximation, novel to this work, which 
is computationally efficient and which has the advantage of being biased in a known way: 
it will always underestimate the amount of shielding, and so we will overestimate the pho- 
todissociation rate of H 2 . The second approach that we have used is our so-called six-ray 
approximation, in which we compute H2 and dust column densities only along the three co- 
ordinate axes of the simulatio n volum e. This approa c h is s imilar to one used previously by 
Nelson fc Langer jl997L Il999h and by lYoshida et al.l (120031 ). but we believe that this is the 



first time that it has been applied in grid-based simulations of this scale. This approximation 
does a much better job of capturing the shielding due to dust, but will tend to overestimate 
the amount of H2 self-shielding once the flow velocities become supersonic. For the simu- 
lations presented in this paper, it proves to be the better choice, but for the simulations 
of turbulent gas presented in paper II, the local shielding approximation proves to be more 
useful. 

Our modifications to ZEUS-MP also include a detailed treatment of the thermal be- 
haviour of the gas. By computing the effects of the most important heating and cooling 
processes active in the ISM, we are able to follow the thermal evolution of the gas with far 
more accuracy than if we relied upon a description in terms of a simple equation of state, 
whether isothermal or polytropic. In dense gas, our accuracy is limited by our inability to 
fully resolve the cooling length of the gas, with the result that we will tend to overproduce 
warm gas. However, the impact of this on the simulations presented in this paper appears 
to be small. 

In § [3j we discussed how the modified code w as tested, and presented the results of a 



test designed to reproduce the IWolfire et al\ (120031 ) analysis of the thermal equilibrium state 



of the neutral atomic ISM. Although we do not reproduce their results exactly - most likely 
due to the fact that the chemical network used in our simulations is very much simpler than 
that adopted in their analysis - we do find qualitatively similar behaviour. Moreover, a 
quantitative comparison of the two sets of results suggests that at most densities, our value 
for the equilibrium gas temperature is accurate to within 50% or better; only for gas at 
densities n ~ 1 cm -3 , which lies in the middle of the thermally unstable regime, do we find 
larger errors, and even here we are accurate to within a factor of two. 

Although our main goal in making these modifications was to study H2 formation in 
turbulent gas (with results that we report on in paper II), we also applied our modified 
hydrodynamical code to the problem of H 2 formation in gravitationally collapsing gas without 
turbulence, starting from initial conditions in which the gas was smoothly distributed and 
initially at rest. We showed that with these initial conditions and an assumed mean density 
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n\ = 100 cm -3 , gravitational collapse occurs on a timescale t co ii > t s . The precise value 
of t coll depends on our choice of shielding approximation (as this significantly affects the 
pressure balance of the gas), and is also sensitive to the size of the simulation volume, the 
strength of the magnetic field, and in some circumstances to the size of the initial density 
inhomogeneities. Our use of periodic boundary conditions probably artificially stabilizes the 
gas to some extent, but even if these were not used, we would still expect collapse to take at 
least one free-fall time (which for our standard initial density rij is approximately 5Myr), and 
in fact probably significantly longer, since the thermal pressure of the gas is never negligible. 

In the simulations that we ran that used the local shielding approximation, we found 
that H2 formation occurred in two phases. During the first phase, which had a timescale of 
5-10 Myr, (xh 2 )m grew slowly until it reached a small limiting value, set by the fact that the 
gas had reached photodissociation equilibrium. Owing to our approximate treatment of H 2 
shielding, this limiting value was resolution dependent (and was smaller in higher resolution 
simulations), but even in our lowest resolution, 64 3 zone simulations, it corresponded to no 
more than about 25% of the gas being molecular. The second, rapid phase of H2 formation 
was triggered by the runaway gravitational collapse of the gas, as the increased density 
boosted the H 2 formation rate while also increasing the amount of shielding against UV 
photodissociation. By the end of the simulations, between 50 and 100% of the gas had 
become molecular, with most of the H 2 having being formed during this collapse phase. 
Therefore, in these simulations, gravitational collapse drives H 2 formation. 

In the simulations that we ran that used the six-ray approximation, the situation was 
somewhat different. In these runs, the amount of shielding in most of the box was con- 
siderably higher, and so the gas generally did not reach photodissociation equilibrium until 
very late times. Therefore, in these runs H 2 formation was reasonably efficient even in the 
absence of dynamical effects. However, we have seen that in these runs the fact that the 
gas in the center is shielded more than the gas at the edges leads to the development of 
a significant temperature gradient and associated pressure gradient. This drives a flow of 
gas towards the center of the box, and the presence of a dynamically significant magnetic 
field causes this flow to be oriented in a direction parallel to the magnetic field lines. The 
resulting pressure-driven compression is ultimately responsible for triggering runaway grav- 
itational collapse. Once this runaway collapse begins, we see a distinct acceleration in the 
rate of H 2 formation, and so even in these simulations it is true to say that gravitational 
collapse plays an important role in driving H 2 formation. These results should be compared 
with the results from simulations of supersonically turbulent gas that we present in paper 
II. They demonstrate that once a sufficiently high mean density is reached, whether through 
gravitational instability or other means, the turbulent compressions alone lead to substantial 
H 2 formation, with gravitational collapse subsequently playing only a minor role. 
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We have also examined the spatial distribution of the H2 produced in our simulations. 
We have shown that by the end of the simulation, most of the H 2 is to be found in dense 
gas, with at least 50% in gas denser than 5000 cm -3 . Moreover, since the initial distributions 
of the gas and the magnetic field are nearly uniform, and since the subsequent dynamical 
evolution proceeds smoothly, the final spatial distribution is highly ordered. Most of the H 2 , 
and, indeed, most of the gas, are found in a thin, dense sheet oriented perpendicularly to 
the direction of the magnetic field. Visually, this distribution looks quite unlike that of the 
molecular gas in real molecular clouds, suggesting that more small-scale power is required in 
either the initial density distribution or velocity distribution (or both) in order to produce 
realistic-looking clouds. 

Finally, we briefly examined the thermal behaviour of gas in our simulations, and showed 
that, as expected, most of the gas is in thermal equilibrium, owing to the short cooling times 
at the simulated gas densities. However, the gas is not isothermal, nor is it describable as a 
simple polytrope: while we can describe the relationship between density and temperature 
with a function of the form 

T oc n 7-1 , (63) 

the polytropic exponent 7 is not constant, but varies with density as shown in Figure [5TJ 
In runs performed using the local shielding approximation, 7 e fj ~ 0.7-0.8 over the range of 
densities covered by our simulations, but in runs that used the six-ray approximation, we 
found a much wider range of 7 c ff. We conclude that while there may be some applications 
for which treating the g simple polytrope with a constant polytropic exponent is a 

reasonable approximation, this approach is probably not valid in the general case. 

Obviously, we do not claim that the results of the study presented in this paper are 
directly applicable to the real ISM: we have investigated a rather simplified dynamical situ- 
ation, which is missing one of the main physical ingredients present in the real ISM, namely 
supersonic turbulence. However, we believe that this study is interesting for the light it 
sheds on the rate at which molecular clouds will form in the absence of turbulence: as we 
have seen, in this case substantial quantities of H 2 are formed only on timescales t > ts, 
consistent with a cloud formation timescale of at least 5-10 Myr. As we shall see in paper 
II, the behaviour of supersonically turbulent gas is very different. 

The authors would like to thank the anonymous referee for suggesting the use of the 
six-ray shielding approximation, and for a number of other comments that have helped us 
to improve the presentation of this paper. They would also like to acknowledge valuable 
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Fig. 12. — Slice in the x-z plane through the density field in 256 3 zone run MS256-RT at a 
time t = 20.6 Myr. The slice is centered on the midpoint of the simulation volume. 
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Fig. 13. — (a) Slice in the x-y plane through the density field in 256 3 zone run MS256 at 
a time t = 28.5 Myr. The periodic boundary conditions employed in the simulation have 
allowed us to shift the image so that the densest region lies at the centre of the figure, (b) 
As (a), but for a slice in the x-z plane. 
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Fig. 14.— (a) Slice in the x-y plane through 256 3 zone run MS256-RT at time t = 17.4 Myr 
showing the spatial variation of the H 2 fraction. The slice is centered on the midpoint of the 
simulation volume, (b) As (a), but for a slice in the x-z plane. 
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Fig. 15. — (a) Slice in the x-y plane through 256 3 zone run MS256 at time t = 28.5 Myr 
showing the spatial variation of the H 2 fraction. The peak value in this slice is 0.77. (b) As 
(a), but for a slice in the x-z plane. 
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Fig. 16. — (a) H2 fraction as a function of the number density of the gas (crosses) for 256 3 
zone run MS256-RT at time t = 17 A Myr. To compute these values, we binned the data 
by number density, using bins of width 0.05 dex, and computed the mean value of xu 2 for 
each bin. The standard deviation in the value of xh 2 in each bin is also indicated (where it 
exceeds the size of the symbol), (b) As (a), but for run MS256 at time t = 28.5 Myr. In this 
figure, we also indicate the mean value of XH 2 ,eq in each bin (dotted line). 
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Fig. 17. — (a) Mass- weighted density PDF, p m (n), in 256 3 zone run MS256-RT at times 
t = 17.4 Myr (solid line) and t = 20.6 Myr (dotted line). Note that gas at n > 5500 cm -3 is 
not properly resolved by the code, and so the gas distribution at t — 20.6 Myr may not be 
quantitatively accurate (although it should be qualitatively correct), (b) As (a), but for run 
MS256 at output times t = 28.5 Myr (solid line) and t = 31.7 Myr (dotted line). 
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Fig. 18. — (a) Cumulative mass distribution of H2 with n in 256 3 zone run MS256-RT at 
times t = 17.4 Myr (solid line) and t = 20.6 Myr (dotted line), (b) As (a), but for run MS256 
at output times t = 28.5 Myr (solid line) and t = 31.7 Myr (dotted line). 
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Fig. 19. — (a) Evolution with time of the minimum and maximum gas temperatures, T min 
and T max , in 256 3 zone runs MS256 (solid lines) and MS256-T100 (dashed lines), which were 
performed with initial temperatures of 1000 K and 100 K respectively, (b) Evolution with 
time of T min and T max in 256 3 zone run MS256-RT. 
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Fig. 20. — (a) Mean gas temperature T plotted as a function of the number density n in 
256 3 zone run MS256-RT at time t = 17.4 Myr. The data were binned as in Figure [TBI above. 
The standard deviation in each bin is also indicated whenever it is larger than the size of 
the symbols used in the plot, (b) As (a), but for run MS256 at time t = 28.5 Myr. 
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Fig. 21.— (a) Value of 7 cff as a function of n in 256 3 zone run MS256-RT at time t = 17.4Myr. 
The data were binned as indicated in the text, (b) As (a), but for run MS256 at time 
t = 28.5 Myr. Note the difference in the vertical scale compared to Figure 21a. 
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Fig. 22. — Evolution with time of the minimum and maximum gas temperatures, T min 
and T max , in 256 3 zone runs MS256 (solid lines) and MS256-CR (dashed lines), which were 
performed with cosmic ray ionization rates ( = 1CT 17 s _1 and ( = 1CT 15 s _1 respectively. 
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Fig. 23. — (a) Evolution of (xu 2 )m with time in a set of 256 3 zone runs in which L was varied 
from 20 pc to 60 pc. The three lines in the upper left represent runs MS256-RT-L20 (solid 
line), MS256-RT (dashed line) and MS256-RT-L60 (dot-dashed line). The five lines in the 
bottom right represent runs MS256-L20 (lower solid line), MS256-L30 (dashed line), MS256 
(dot-dashed line), MS256-L50 (dotted line) and MS256-L60 (upper solid line), (b) As (a), 
but for the maximum gas density, p max . 
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Fig. 24. — Evolution of (xh 2 )m with time in a set of 256 3 zone runs in which the maximum 
amplitude of the initial density perturbations, 5, was varied. Results are plotted for runs 
MS256-RT and MS256-RT-dl00 (solid and dashed lines in the upper left) as well as for runs 
MS256 (solid line), MS256-dlO (dashed line), MS256-d50 (dot-dashed line) and MS256-dl00 
(dotted line). The first pair of runs have 5 = 0.05 and 1.0 respectively; the latter four have 
5 = 0.05, 0.1, 0.5 and 1.0 respectively. 
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Fig. 25. — (a) Evolution of (xu 2 )m with time in a set of 256 3 zone runs in which the strength 
of the initial magnetic field was varied. Results are plotted for seven runs. Three of these 
runs used the six-ray shielding approximation: runs HS256-RT (upper solid line), MS256-RT 
(upper dashed line) and MS256-RT-Bx4 (upper dot-dashed line), which had initial magnetic 
field strengths B\ = 0.0, 5.85 and 23.4 fiG respectively. The other four runs used the local 
shielding approximation: runs HS256 (lower solid line), MS256 (lower dashed line), MS256- 
Bx2 (lower dot-dashed line) and MS256-Bx4 (lower dotted line), which had initial magnetic 
field strengths B\ = 0.0, 5.85, 11.7 and 23.4 /j,G respectively. In the magnetized runs, there 
is so little sensitivity to B\ that the lines are hard to distinguish from each other in the plot. 
On the other hand, in the B x — runs we see a significant difference in behaviour, (b) As 
(a), but for the maximum gas density p max - 
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Fig. 26. — (a) Evolution of (:th 2 )m with time in a set of 256 3 runs in which the initial gas 
density was varied. Results are plotted for runs MS256-RT (lefthand solid line), MS256- 
RT-nlO (lefthand dashed line), MS256 (righthand solid line), MS256-n30 (righthand dashed 
line) and MS256-nlO (dot-dashed line). Note that L and were varied in these runs so 
as to keep the number of Jeans masses in the simulation volume and the mass-to-flux ratio 
approximately constant, (b) As (a), but for p max - 
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Table 1. The set of chemical reactions that make up our model of non-equilibrium 

hydrogen chemistry. 

Reaction Reference 

1. H + H + grain — >■ H 2 + grain Hollenbach fc McKee (1979) 

2. H 2 + H -> 3H Mac Low fc Shull (1986) (low density), 

Lepp fc Shull (1983) (high density) 

3. H 2 + H 2 -> 2H + H 2 Martin. Keogh & Mandv (1998) (low density), 

Shapiro fc Kang (1987) (high density) 

4. H 2 + 7 -> 2H See § I2T2TT1 

5. H + c.r. -> H+ + e See § O 

6. H + e -> H+ + 2e Abel et al. (1997) 

7. H+ + e -> H + 7 Ferland of. (1992) 

8. H + + e + grain — > H + grain Weingartner fc Draine (2001) 
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Table 2. Processes included in our thermal model. 



Process 



References 



Cooling: 

C II fine structure lines 



O I fine structure lines 



Si II fine structure lines 



H2 rovibrational lines 
Gas-grain energy transfer 1 
Recombination on grains 
Atomic resonance lines 
H collisional ionization 
H2 collisional dissociation 



Atomic data - Silva fc Viegas (2002) 

Collisional rates (H 2 ) - Flower fc Launav ( 19771 

Collisional rates (H, T < 2000 K) - Hollenbach fc McKee (1989) 

Collisional rates (H, T > 2000 K) - Keenan et al. (1986) 

Collisional rates (e~) - Wilson fc Bell (2002) 

Atomic data - Silva fc Viegas (2002) 

Collisional rates (H, H 2 ) - Flower, priv. comm. 

Collisional rates (e~) - Be ll. Berrington fc Thoma s (1998) 

Collisional rates (H+) - Peauignot (1990 , 1996) 

Atomic data - Silva fc Viegas (2002) 

Collisional rates (H) - Roueff (1990) 

Collisional rates (e~) - Dufton fc Kingston (1991) 

Le Bourlot. Pineau des Forets fc Flower (1999) 

Hollenbach fc McKee (1989) 

Wolfire et al. (2003) 

Sutherland fc Dopita (1993) 

Abel et al. (1997) 

See Tabled] 



Heating: 

Photoelectric effect 

H 2 photodissociation 

UV pumping of H 2 

H 2 formation on dust grains 

Cosmic ray ionization 



Bakes & Tielens (1994) : Wolfire et al. (2003) 

Black fc Dakarno (1977) 

Burton. Hollenbach fc Tielens (1990) 

Hollenbach fc McKee (1989) 

Goldsmith fc Langer (1978) 



Note. - 1: If T gas < T grain , the net flow of energy is from the grains to the gas, leading to 
heating instead of cooling. 
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Table 3. The number density at which the Truelove criterion is violated, n max , computed 
for gas in thermal and chemical equilibrium for various different box sizes and numerical 

resolutions. 



Resolution 


Box size (pc) 


Zone size (pc) 


n max (cm 3 ) 




20 


0.31 


2.0 x 10 3 


64 3 


40 


0.62 


8.0 x 10 2 




60 


0.94 


4.9 x 10 2 




20 


0.16 


5.4 x 10 3 


128 3 


40 


0.31 


2.0 x 10 3 




60 


0.47 


1.2 x 10 3 




20 


0.078 


1.5 x 10 4 


256 3 


40 


0.16 


5.4 x 10 3 




60 


0.23 


3.0 x 10 3 




20 


0.039 


5.6 x 10 4 


512 3 


40 


0.078 


1.5 x 10 4 




60 


0.12 


8.2 x 10 3 
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Table 4. Input parameters used for our runs 



Run 


L he) 


5 rii (cm 3 ) 


Ti (10 


Bi (uG) 


Notes 


MS64 


40 


0. 


.05 


100 


1000 


5.85 




MS128 


40 





,05 


100 


1000 


5.85 




MS256 


40 





,05 


100 


1000 


5.85 




MS64-RT 


40 


o 


,05 


100 


1000 


5.85 




MS128-RT 


40 


0. 


.05 


100 


1000 


5.85 




MS256-RT 


40 


0. 


.05 


100 


1000 


5.85 




MS64-ng 


40 





.05 


100 


1000 


5.85 


1 


MS128-ng 


40 


0. 


.05 


100 


1000 


5.85 


1 


MS256-ng 


40 





.05 


100 


1000 


5.85 


1 


MS256-RT-ng 


40 


0. 


.05 


100 


1000 


5.85 


1 


MS64-nr 


40 


o 


,05 


100 


1000 


5.85 


2 


MS128-nr 


40 


0. 


.05 


100 


1000 


5.85 


2 


MS256-nr 


40 


o 


.05 


100 


1000 


5.85 


2 


MS256-T100 


40 


0. 


.05 


100 


100 


5.85 




MS256-CR 


40 





.05 


100 


1000 


5.85 


3 


MS256-L20 


20 





.05 


100 


1000 


5.85 




MS256-RT-L20 


20 


0. 


.05 


100 


1000 


5.85 




MS256-L30 


30 





.05 


100 


1000 


5.85 




MS256-L50 


50 


0. 


.05 


100 


1000 


5.85 




MS256-L60 


60 


n 
u 


,uo 


1 nn 
ruu 


1000 


5.85 




MS256-RT-L60 


60 


0. 


.05 


100 


1000 


5.85 




MS256-Bx2 


40 





,05 


100 


1000 


11.7 




MS256-Bx4 


40 


0. 


.05 


100 


1000 


23.4 




MS256-RT-Bx4 


40 





,05 


100 


1000 


23.4 




HS256 


40 





.05 


100 


1000 


0.0 




HS256-RT 


40 


0. 


.05 


100 


1000 


0.0 




MS256-dlO 


40 


0. 


.10 


100 


1000 


5.85 




MS256-d50 


40 


0. 


.50 


100 


1000 


5.85 




MS256-dl00 


40 


1 


,00 


100 


1000 


5.85 




MS256-RT-dl00 


40 


1 


.00 


100 


1000 


5.85 




MS256-nlO 


85 





,05 


10 


1000 


1.24 
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Table 4 — Continued 



Run 


L (pc) 


5 


rii (cm 3 ) 


T, (K) 


Bi (/iG) Notes 


MS256-RT-nlO 


85 


0.05 


10 


1000 


1.24 


MS256-n30 


60 


0.05 


30 


1000 


2.63 



Note. — 1: runs with self-gravity disabled; 2: runs with x — 0.0; 3: run 
with a higher cosmic ray ionization rate, ( = 10^ 15 s _1 . 
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Table 5. t rcs , £f, and associated values of (xh 2 )m for all runs in Table HI 



Run t res (Myr) (xj^M^res) U (Myr) (x-B. a )u{tf) 



MS64 


21.3 


0.30 


27.3 


0.89 


MS128 


26.5 


0.34 


31.1 


0.93 


MS256 


29.5 


0.21 


31.7 


0.50 


MS64-RT 


15.0 


0.65 


22.1 


0.99 


MS128-RT 


16.9 


0.74 


20.9 


0.99 


MS256-RT 


18.4 


0.89 


20.6 


0.99 


MS64-ng 


— 


— 


31.7 


0.26 


MS128-ng 


— 


— 


31.7 


0.13 


MS256-ng 


— 


— 


31.7 


6.4 x 10 


MS256-RT-ng 


— 


— 


19.9 


0.73 


MS64-nr 


18.3 


0.70 


29.6 


0.99 


MS128-nr 


22.1 


0.76 


28.9 


0.98 


MS256-nr 


26.6 


0.83 


31.7 


0.96 


MS256-T100 


28.9 


0.23 


31.7 


0.62 


MS256-CR 


30.0 


0.42 


31.7 


0.76 


MS256-L20 


— 


— 


31.7 


7.4 x 10 


MS256-RT-L20 


17.2 


0.96 


20.0 


0.99 


MS256-L30 


31.3 


0.42 


31.7 


0.48 


MS256-L50 


27.7 


0.20 


30.8 


0.71 


MS256-L60 


26.7 


0.20 


29.9 


0.70 


MS256-RT-L60 


17.8 


0.76 


22.1 


0.98 


MS256-Bx2 


29.8 


0.25 


31.7 


0.50 


MS256-Bx4 


29.2 


0.27 


31.7 


0.62 


MS256-RT-Bx4 


18.7 


0.92 


19.0 


0.96 


HS256 


23.9 


0.10 


25.2 


0.21 


HS256-RT 


13.1 


0.60 


18.2 


0.97 


MS256-dlO 


27.1 


0.22 


30.5 


0.68 


MS256-d50 


22.9 


0.20 


25.8 


0.64 


MS256-dl00 


21.0 


0.19 


23.7 


0.59 


MS256-RT-dl00 


17.1 


0.80 


20.2 


0.99 


MS256-nlO 






95.1 


8.7 x 10 
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Table 5 — Continued 



Run 


t rcs (Myr) 


(sh 2 )m(*i«j) U (Myr) (x H2 ) M (tf) 


MS256-RT-nlO 




53.0 0.91 


MS256-n30 


58.1 


0.26 64.1 0.70 



Note. - t res is the time at which the Truelove criterion is first 
violated during the course of the run; when no value is given, this 
indicates that the criterion was never violated, tf is the time at which 
the simulation was stopped. 



